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Abstract 

This  study  presents  a  method  of  moments  solution  for  the  surface 
currents  and  charge  distributions  on  an  aperture-fed,  stacked  patch 
antenna.  To  make  the  solution  independent  of  the  antenna  excitation 
technique,  an  aperture  magnetic  current  distribution  is  assumed,  and  the 
spatial  Green's  functions  of  the  antenna  are  used  to  calculate  the 
tangential  fields  on  the  antenna  patches  and  coefficients  for  a  finite 

series  expression  of  the  surface  currents  and  charge  densities.  The 

Green's  functions  of  the  antenna  are  modeled  with  polynomials  that  are  a 

function  of  the  radial  distance  separating  observer  and  source  loca¬ 
tions.  The  method  of  moments  solution  is  implemented  using  FORTRAN  77 
and  requires  no  external  software  support.  The  numerical  techniques, 
algorithms,  and  complete  software  functional  descriptions  are  also 
included. 
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A  METHOD  OF  MOMENTS  SOLUTION  FOR  THE  ELECTRIC 


CURRENTS  ON  AN  APERTURE-FED,  STACKED  PATCH  MICROSTRIP  ANTENNA 


I.  Introduction 


This  chapter  presents  a  brief  overview  of  microstrip  antennas  and 
the  motivation  for  the  analysis  of  the  aperture-fed,  stacked  patch 
microstrip  antenna.  The  exact  goals  of  this  thesis,  along  with  the  as¬ 
sumptions  made  during  the  research,  are  discussed.  A  brief  review  of 
current  aperture-fed,  stacked  patch  microstrip  antenna  analysis  efforts, 
and  an  introduction  to  the  method  of  moments  is  also  included.  In  con¬ 
clusion,  an  overview  of  the  remaining  chapters  of  the  thesis  is  given. 


1.1  Background  on  Microstrip  Antennas 

In  its  simplest  form,  a  microstrip  antenna  consists  of  two  paral¬ 
lel,  conducting  layers  separated  by  a  dielectric  substrate.  Typically, 
the  upper  layer  is  the  radiating  element  and  the  lower  serves  as  a 


ground  plane.  Input  to  the  antenna  can  be  supplied  in  a  number  of  ways. 
Two  of  the  most  versatile  feed  techniques  are  the  microstrip  feedline 
and  the  coaxial  feed  (see  Figure  1-1).  Radiation  from  the  antenna  is 


Microstrip 

Feedline 


Radiating  Patches 


Ground  Plane 

4 

Coaxial  Feed 


(a) 


(b) 


Figure  1-1  Microstrip  Antenna  with  (a)  Microstrip  Feed  and  (b)  Coaxial 
Feed 


produced  by  fringing  of  the  electric  field  along  the  open-circuited 
edges  of  the  patch.  The  radiating  patch  of  the  antenna  is  typically  A/2 
long,  thus  the  normal  components  of  the  antenna's  radiation  tend  to 
cancel  in  the  far-field,  while  the  components  parallel  to  the  antenna 
ground  plane  are  in  phase.  Therefore,  the  radiated  fields  are  at  a 
maximum  in  the  direction  normal  to  the  antenna's  surface  [1:5).  The 
radiating  patch  may  be  any  shape,  but  rectangular  or  circular  elements 
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are  commonly  used  because  of  the  relative  simplicity  of  analyzing  their 
radiation  characteristics  [2:488]. 

G.  A.  Deschamps  first  proposed  the  concept  of  a  microstrip  antenna 
in  1953.  However,  20  years  passed  before  photo-etching  techniques  and 
material  properties  were  developed  to  the  point  where  practical  micro¬ 
strip  antennas  could  be  built  [1:11.  A  series  of  successful,  indepen¬ 
dent  experiments  in  the  early  1970 's  demonstrated  the  feasibility  of 
microstrip  antennas.  By  1978,  the  microstrip  antenna  was  widely  known 
and  used  in  a  variety  of  communications  systems  [6:2].  Since  the 
1970' s,  extensive  research  and  development  exploiting  the  numerous  ad¬ 
vantages  of  the  microstrip  antenna,  both  as  a  single  radiation  source 
and  as  an  array  element,  have  established  it  as  a  separate  entity  within 
the  broad  field  of  microwave  antennas  [1:1]. 


1.2  Advantages  and  Disadvantages  of  Microstrip  Antennas 

Microstrip  antennas  are  well  suited  to  the  current  emphasis  on 
miniaturization  and  increased  reliability.  The  antennas  have  been 
successfully  adapted  for  use  in  systems  ranging  from  biomedical  diagno¬ 
sis  to  satellite  and  aircraft  communications  [20:424], 

The  widespread  use  of  microstrip  antennas  can  be  attributed  to  the 
numerous  advantages  they  offer  when  compared  to  more  traditional  antenna 
types.  Some  of  these  advantages  include  low  volume,  weight,  and  fabri¬ 
cation  cost,  mechanical  ruggedness,  and  compatibility  with  modular 
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designs,  since  solid  state  devices  can  be  manufactured  simultaneously 
and  on  the  same  circuit  board  as  the  antenna  [23:2]. 

Some  of  the  disadvantages  associated  with  microstrip  antennas 
include  narrow  bandwidth,  antenna  gains  much  less  than  unity,  poor 
isolation  between  feed  and  radiating  elements,  and  low  power  handling 
capability  [1:2]. 

Narrow  bandwidth  is  the  most  restrictive  of  these  shortcomings. 

The  bandwidth  of  a  simple  microstrip  antenna  is  typically  less  than  five 
percent  of  the  resonant  frequency  [21:145].  The  bandwidth  of  a  micro¬ 
strip  antenna  can  be  improved  by  reducing  the  permittivity  or  increasing 
the  thickness  of  the  artenna  substrate  [1:62].  However,  using  a  low 
permittivity  substrate  degrades  the  performance  of  other  solid  state 
components  etched  on  the  same  substrate  as  the  antenna.  Increasing  sub¬ 
strate  thickness  detracts  from  the  microstrip  antenna's  low  profile,  and 
causes  difficulties  in  matching  the  antenna  input  impedance  with  that  of 
the  feed  network  [23:3]. 


1.3  The  Aperture-Fed ,  Stacked  Patch  Microstrip  Antenna 

Several  practical  ways  to  overcome  bandwidth  limitations  without 
compromising  the  advantages  of  the  microstrip  antenna  have  been  tried 
and  proven.  The  aperture-fed,  stacked  patch  antenna  (see  Figure  1-1)  is 
effort  to  improve  bandwidth  by  employing  two  of  these  techniques,  the 
aperture-fed  and  stacked  patch  configurations,  in  a  single  design. 


Antenna 

By  placing  the  antenna  patches  on  a  substrate  separated  from  other 
active  circuit  components,  the  aperture-fed  configuration  allows 
selection  of  substrate  permittivity  and  thickness  for  maximum  antenna 
performance  without  affecting  the  operation  of  other  circuit  elements. 
Additionally,  placing  the  ground  plane  between  the  feed  network  and  the 
radiating  elements  of  the  antenna  overcomes  the  problem  of  poor  isola¬ 
tion  between  spurious  feed  network  and  antenna  radiation  [31:936]. 
Because  there  is  no  physical  connection  between  the  feed  network  and  the 
antenna,  the  aperture-fed  configuration  is  simple  to  manufacture  and 
consequently,  very  appealing  for  use  in  antenna  arrays,  especially  those 
containing  a  large  number  of  elements. 


1-5 


The  use  of  the  stacked  patch  configuration  provides  further 
bandwidth  enhancement  [23:3].  Proper  selection  of  the  additional 
dielectric  layer's  permittivity  and  thickness,  as  well  as  the  shapes, 
sizes,  and  relative  locations  of  the  feed  and  radiating  patches,  can 
provide  the  additional  benefits  of  enhanced  gain  or  dual  frequency 
operation  [5:609]. 


1.4  Research  Description 

The  aperture-fed,  stacked  patch  microstrip  antenna  presents  a 
complex  analysis  problem.  One  analysis  has  been  completed  by  Nazar  [23] 
using  the  mixed  potential  integral  equation  (MPIE)  approach.  Nazar's 
work  provides  the  integral  equations  that  can  be  used  to  calculate  the 
currents  and  potentials  on  an  aperture-fed,  stacked  patch  microstrip 
antenna.  However,  the  integral  equations  cannot  be  solved  analytically; 
thus,  a  numerical  solution  is  required. 

1.4.1  Problem  Statement.  Develop  the  FORTRAN  software  necessary 
to  calculate  the  electric  surface  currents  and  charge  densities  on  the 
elements  of  an  aperture-fed,  stacked  patch  microstrip  antenna. 

1.4.2  Research  Objectives  and  Scope.  This  effort  uses  Nazar's 
analysis  as  a  basis  for  developing  a  Method  of  Moments  (MoM)  solution 
for  the  surface  currents  and  charge  densities  on  an  aperture-fed, 
stacked  patch  microstrip  antenna.  With  this  solution,  the  radiated 
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fields  and  impedances  characteristics  of  the  antenna  design  can  be 
calculated. 

The  product  of  this  research  is  a  software  package  that  completes 
the  MoM  solution  for  a  user  specified  aperture-fed,  stacked  patch 
microstrip  antenna.  The  package  is  written  in  FORTRAN  77  and  requires 
no  external  software  to  complete  the  solution.  Results  obtained  with 
the  package  have  been  verified  as  much  as  possible  using  data  available 
in  current  literature  [18].  Descriptions  explaining  the  use  of  the 
software  and  the  design,  coding,  and  logic  of  the  solution  are  also  pro¬ 
vided. 

1.4.3  Research  Assumptions .  Nazar's  analysis  uses  Green's  func¬ 
tions  to  characterize  the  interactions  and  radiation  mechanisms  of  the 
aperture-fed,  stacked  patch  antenna.  By  definition,  the  Green’s 
functions  are  expressions  for  the  potentials  created  by  a  unit  electric 
or  Hertz  dipole  located  within  one  of  the  antenna's  surfaces  or  ground 
plane.  By  expressing  any  arbitrary  signal  applied  to  the  antenna  as  a 
spectrum  of  unit  sources,  the  Green's  functions  can  be  used  to  calculate 
the  resulting  fields  in  the  antenna  [19:165],  With  the  fields  known,  it 
is  then  possible  to  calculate  the  radiated  fields. 

Nazar  expresses  the  Green's  functions  as  integral  equations  that 
are  derived  under  the  restrictions  of  each  dielectric  layer  having 
finite  thickness  while  being  homogeneous,  isotropic,  and  lossless. 
Additionally,  the  patch  conductors  and  ground  plane  are  infinitely  thin, 
perfect  conductors. 
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The  software  developed  for  this  thesis  does  not  explicitly 
evaluate  the  integral  expressions  for  the  aperture-fed,  stacked  patch 
microstrip  antenna's  Green's  functions.  Instead,  the  Green's  functions 
are  represented  by  appropriately  chosen  polynomials  in  accordance  with 
the  work  of  Erwert  [8]. 


1.5  Summary  of  Current  Knowledge 

As  previously  discussed,  Nazar  has  derived  the  Green's  Functions 
necessary  to  model  the  aperture-fed,  stacked  patch  microstrip  antenna. 
His  work  combines  the  aperture-feed  analysis  of  Sullivan  and  Schaubert 
[29]  with  the  MPIE  approach  of  Mosig  and  Gardiol  [20]  to  derive  the 
spatial  Green's  functions  for  the  antenna.  The  MPIE  approach  is  attrac¬ 
tive  because  it  fully  characterizes  the  near-fields  of  the  antenna  in¬ 
cluding  surface  waves  [ 23 : ii ] .  The  Green's  functions  for  the  antenna 
are  expressed  in  the  form  of  Sommerfeld  integrals  that  can  only  be 
evaluated  numerically.  However,  the  numerical  stability  of  the  Green’s 
functions  allow  the  use  of  approximations  for  their  values  without 
unacceptable  loss  of  accuracy  [ 19 : 429 ; 21 : 224-226 ] 

The  method  of  moments  is  a  numerical  technique  that  has  been 
successfully  used  to  solve  a  wide  variety  of  electromagnetic  radiation 
and  scattering  problems.  The  MoM  approach  provides  a  well  understood 
means  for  converting  a  series  of  complex,  integral  equations  into  matrix 
equations  that  are  easily  solved  on  a  digital  computer.  Furthermore, 
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the  MoM  technique  has  been  used  in  microstrip  antenna  radiation  solu¬ 
tions  by  numerous  authors  [16;22;24],  Sullivan  and  Schaubert  [29]  have 
successfully  applied  the  method  of  moments  in  the  analysis  of  an 
aperture-fed  antenna,  and  Barlatey  et  al.  [5]  have  employed  the  tech¬ 
nique  to  obtain  the  MPIE  solution  of  the  stacked  patch  configuration.' 


1.6  Materials  and  Equipment  Needed 

Theoretically,  the  MoM  software  package  can  be  executed  on  any 
computer  system  equipped  with  a  FORTRAN  77  compiler.  However,  due  to 
the  large  number  of  numerical  integrations  required  by  the  MoM  solution, 
a  high  speed  computer  is  best  suited  for  the  task.  In  this  study,  the 
calculations  were  accomplished  on  the  Air  Force  Institute  of  Technology, 
School  of  Engineering  Digital  Equipment  Corporation  VAX  8550  computer. 


1.7  Thesis  Organization 

The  remainder  of  this  document  is  organized  as  follows.  Chapter 
II  presents  a  survey  of  literature  focusing  on  current  research  efforts 
involving  MPIE  analysis  of  aperture-fed  microstrip  antennas,  stacked 
patch  microstrip  antennas,  or  aperture-fed,  stacked  patch  microstrip 
antennas  and  the  use  of  the  method  of  moments  in  obtaining  solutions  to 
microstrip  antennas.  Chapter  III  presents  the  equations  needed  to  model 


\ 
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the  antenna  and  the  specifics  of  the  MoM  solution.  Chapter  IV  discusses 
the  numerical  techniques  needed  to  complete  the  MoM  solution  and  a 
detailed  description  of  the  logic,  algorithms,  and  FORTRAN  code. 

Chapter  V  discusses  the  predictions  made  with  the  software  and  their 
accuracy  as  compared  to  data  reported  in  the  literature.  Finally,  the 
conclusions  drawn  from  this  research  and  recommendations  for  further 
work  are  presented  in  Chapter  VI. 
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II.  Current  Efforts  -  MLcrostrip  Antennas 


This  chapter  presents  a  review  of  current  microstrip  analysis 
efforts  relevant  to  the  research  of  this  thesis.  To  begin,  a  detailed 
discussion  of  aperture-fed  and  stacked  patch  microstrip  antennas  is 
presented  with  emphasis  on  the  benefits  of  the  configurations.  Next, 
research  on  aperture-fed,  stacked  patch  antennas  is  reviewed,  and 
finally,  applying  the  method  of  moments  to  solutions  for  microstrip 
antennas  is  discussed. 


2.1  Aperture  Fed  Microstrip  Antennas 

The  bandwidth  performance  of  a  microstrip  antenna  can  be  improved 
by  either  increasing  the  thickness  or  decreasing  the  relative  permit¬ 
tivity  of  the  dielectric  substrate  separating  the  conductors  of  the 
antenna  [1:62].  However,  these  changes  result  in  loss,  radiation,  and 
input  impedance  mismatch  problems  which  would  typically  offset  any 
bandwidth  gains  [23:3].  One  way  to  take  advantage  of  increased  sub- 
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strate  thickness  while  avoiding  these  problems  is  the  aperture-fed 
microstrip  antenna  which  was  first  proposed  by  D.  M.  Pozar  [26]. 

Pozar’s  configuration  places  a  ground  plane  between  a  microstrip 
feedline  and  a  microstrip  antenna  patch  (see  Figure  2-1).  The  signal 


(a)  (b) 


Figure  2-1  Side  View  (a)  and  Top  View  (b)  of  Pozar' s  Aperture-Fed 
Microstrip  Antenna 


applied  to  the  microstrip  feedline  is  electromagnetically  coupled  to  the 

radiating  patch,  thus  the  aperture-fed  configuration  does  not  require  a 

direct  connection  between  the  feedline  and  antenna.  Pozar  notes  three 

advantages  to  this  configuration: 

(i)  The  configuration  is  well  suited  for  monolithic  phased 
arrays,  where  active  devices  can  be  integrated  on,  for 
example,  a  gallium  arsenide  substrate  with  the  feed  network 
and  the  radiating  elements  can  be  located  on  an  adjacent 
(low  dielectric  constant)  substrate,  and  coupled  to  the  feed 
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network  through  apertures  in  the  ground  plane  separating  the 
two  substrates.  The  use  of  two  substrates  thus  avoids  the 
deleterious  effect  of  a  high-dielectric-constant  on  the 
bandwidth  and  scan  performance  of  a  printed  antenna  array. 

(ii)  No  radiation  from  the  feed  network  can  interfere  with 
the  main  radiation  pattern,  since  a  ground  plane  separates 
the  two  mechanisms. 

(iii)  No  direct  connection  is  made  to  the  antenna  elements, 
so  problems  such  as  large  probe  self-reactances  or  wide 
microstrip  line  (relative  to  patch  size),  which  are  critical 
at  millimeter-wave  frequencies,  are  avoided.  [26:49] 

Pozar  uses  the  cavity  model  of  a  patch  antenna  and  small  hole 

coupling  theory  to  analyze  the  aperture  coupling  and  quantify  the 

relationship  between  aperture  placement  and  the  coupling  of  energy  to 

the  radiating  patch.  Pozar  found  an  aperture  with  a  long  y-dimension 

(see  Figure  2-1)  and  a  narrow  x-dimension  provides  the  strongest 

coupling  between  the  feedline  and  radiating  patch  [25:50],  Also, 

varying  the  length  of  the  feedline  extending  beyond  the  aperture  and  the 

aperture  size  provides  two  means  of  tuning  the  antenna  [25:50].  Pozar 

suggests  the  use  of  as  small  an  aperture  as  possible,  while  adjusting 

the  feedline  length  to  match  the  antenna  input  impedance  with  that  of 

the  feed  circuit.  While  no  bandwidth  performance  is  reported,  Pozar 

does  note  the  design  appears  to  be  quite  robust  in  terms  of  flexibility 

and  ease  of  fabrication  [25:49].  Measurements  of  a  prototype  antenna 

employing  the  aperture  feed  show  the  principal  plane  radiation  patterns 

are  basically  the  same  as  other  microstrip  antennas  of  the  same  geometry 

[25:50] . 

Sullivan  and  Schaubert  [29]  have  completed  a  more  detailed  math¬ 
ematical  analysis  of  an  aperture-fed  microstrip  antenna  with  a  single 


2-3 


patch.  Their  analysis  uses  coupled  integral  equations  and  a  Galerkin 
moment  method  to  solve  for  the  feedline  electrical  surface  currents,  the 
radiating  patch  electrical  surface  currents,  and  the  equivalent  magnetic 
surface  currents  on  the  ground  plane.  As  part  of  this  work,  Sullivan 
and  Schaubert  used  their  model  to  investigate  the  effect  of  increased 
feedline  thickness  and  relative  permittivity.  They  report  increased 
relative  permittivity  of  the  feedline  substrate  improved  coupling,  but 
left  the  resonant  frequency  of  the  antenna  unchanged  [29:982].  As 
feedline  substrate  thickness  was  increased,  aperture  coupling  decreased 
while  the  resonant  frequency  of  the  antenna  remained  unchanged  [29:982], 
These  results  demonstrate  the  wide  range  of  design  flexibility  offered 
by  the  aperture-fed  antenna.  Additionally,  the  analysis  provides  a 
means  to  predict  antenna  input  impedance  and  resonant  frequency  based  on 
aperture  length  and  the  thickness  and  dielectric  constant  of  the  feed¬ 
line  and  antenna  substrates. 


2.2  Stacked  Patch  Microstrip  Antennas 

An  example  of  a  coaxially-fed,  stacked  patch  microstrip  antenna  is 
shown  in  Figure  2-2.  The  lower  patch  of  the  antenna  is  referred  to  as 
the  driven  or  feed  patch.  The  upper  patch  is  the  radiating  element  and 
is  excited  by  energy  electromagnetically  coupled  from  the  driven  patch. 
Because  of  this  coupling  mechanism,  stacked  patch  antennas  are  also 
referred  to  as  electromagnetically  coupled  patch  (EMCP)  antennas. 


Figure  2-2  Coaxially-Fed,  Stacked  Patch  Microstrip  Antenna 


Oltman  and  Huebner  [25]  first  introduced  the  concept  of  electro- 
magnetically  coupled  microstrip  radiators  with  a  study  of  microstrip 
dipoles  stacked  near  a  microstrip  feedline. 

The  study  cites  four  advantages  offered  by  the  stacked  configura¬ 
tion.  Because  of  increased  ground  plane-radiator  spacing,  the  stacked 
geometry  provides  greater  bandwidth  performance  and  increased  radiation 
efficiency  [25:151].  Experimental  data  reported  in  the  study  show 
bandwidth  values  between  2.5  and  5.5  percent  for  VSWR  of  1.92  or  higher 
[25:154] . 

A  third  advantage  of  the  stacked  geometry  is  reduced  radiation 
from  feedline  bends,  junctions,  and  other  discontinuities.  In  the 
stacked  geometry,  the  radiating  element  can  be  placed  farther  away  from 
the  feed  circuitry,  while  the  feed  circuitry  can  be  simultaneously 
located  close  to  the  ground  plane.  Since  feedline  radiation  is  directly 
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proportional  to  the  separation  between  the  feedline  and  ground  plane, 
the  feedline  radiation  is  reduced  [25:152]. 

The  fourth  advantage  of  the  stacked  geometry  is  the  ease  in 
matching  the  electromagnetically  coupled  radiator  to  its  feedline.  For 
a  basic  microstrip  radiator,  the  radiator  configuration  that  provides  a 
perfect  match  is  unique.  Since  the  stacked  geometry  employs  electric 
and  magnetic  field  coupling,  Oltman  and  Huebner  describe  the  perfect 
match  configuration  as  a  near  ellipsoidal  surface,  .'hereforc  the  perfect 
match  point  is  easier  to  achieve  [25:152]. 

Sabban  [28]  expands  on  the  stacked  element  concept  by  replacing 
the  microstrip  dipoles  with  patches.  In  Sabban' s  design,  the  lower 
patch  is  sized  to  resonate  at  center  frequency,  and  the  upper  patch  is 
slightly  larger  than  the  lower.  A  low  permittivity  material,  such  as 
air,  is  used  to  separate  the  patches.  Several  antennas  were  constructed 
and  tested  as  single  radiators  and  as  array  elements  over  the  S  to  Ka 
frequency  bands.  Sabban  reports  adjustments  to  patch  separation  enable 
impedance  matching  between  the  antenna  and  the  feed  circuit  over  a  wide 
frequency  range  [28:64].  Depending  on  the  antenna  geometry,  bandwidths 
of  9  to  15  percent  with  VSWR  less  than  2:1  were  obtained  [28:64], 


2.3  Aperture  Fed,  Stacked  Patch  Antennas 


Tsao  et  al  [31]  combine  the  aperture  feed  and  stacked  patch 
configurations  to  achieve  even  larger  antenna  bandwidths.  Their  purpose 
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in  using  the  aperture-fed  antenna  is  two  fold.  First,  the  aperture-fed 
antenna  is  simpler  to  manufacture  since  no  probe  must  be  soldered  to  the 
antenna.  Second,  the  ground  plane  separating  the  feed  network  and  the 
antenna  prevents  interference  between  spurious  feed  network  radiation 
and  antenna  radiation  [31:936].  Proper  selection  of  the  thicknesses  and 
permittivities  of  dielectric  layers  separating  the  antenna  elements, 
employed  in  conjunction  with  appropriately  selected  geometries  for  the 
coupling  aperture  and  microstrip  patches,  provides  a  means  to  meet 
various  antenna  bandwidths  requirements  [31:936],  Depending  on  desired 
polarization  characteristics,  bandwidths  of  19.2  and  23  percent  were 
obtained  for  VSWR  less  than  2  [31:936], 

In  [23]  ,  Nazar  combines  the  MPIE  and  stratified  media  techniques 
to  perform  a  detailed  analysis  of  the  aperture-fed,  stacked  patch 
antenna.  Using  the  vector  and  scalar  potentials  for  infinitesimal 
sources  embedded  in  stratified  media,  Nazar  derives  a  set  of  Green's 
functions  in  the  form  of  Sommerfeld  integrals  expressed  in  terms  of 
antenna  spatial  variables  [23:13-39].  Nazar  describes  the  characteris¬ 
tics  of  the  Green's  functions'  integrands,  discusses  integrand  problem 
areas,  and  provides  FORTRAN  software  for  numerical  evaluation  of  the 
integrals  [23:49-84],  Also,  a  technique  for  using  the  Green's  func¬ 
tions  in  a  method  of  moments  solution  for  the  antenna  currents  is 
described  [23:39-48]. 
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2. 4  Microstrip  Antenna  Solutions  Using  the  Method  of  Moments 

As  previously  mentioned,  the  method  of  moments  is  a  numerical 
technique  widely  used  to  solve  electromagnetic  problems.  The  MoM 
provides  a  well  understood  and  easily  implemented  means  of  transforming 
linear  functional  equations,  i.e.  integral,  differential,  or  integro- 
differential  equations,  into  a  matrix  equation  easily  solved  on  a 
computer.  The  equation  being  transformed  typically  describes  an 
appropriately  selected  boundary  condition  that  must  be  satisfied.  In 
the  case  of  a  microstrip  antenna,  the  boundary  condition  is  the  require¬ 
ment  of  all  tangential  electric  field  components  summing  to  zero.  While 
the  MoM  technique  is  a  unique  process,  the  methods  by  which  each  step  of 
the  process  is  completed  are  diverse  and  bear  careful  consideration. 

This  review  concentrates  on  past  efforts  that  provide  insight  into 
choosing  the  features  which  provide  the  best  compromise  between  accuracy 
and  computational  feasibility  of  the  MoM  solution. 

As  employed  in  microstrip  antenna  problems,  the  MoM  breaks  the 
scattering  or  radiating  body  into  a  collection  of  smaller  bodies, 
typically  referred  to  as  cells.  A  simple  approximation  for  the  surface 
currents  on  each  cell,  referred  to  as  a  basis  function,  is  proposed. 
Therefore,  the  surface  currents  of  the  complete  body  is  the  finite  sum 
of  appropriately  weighted  basis  functions  of  each  cell.  The  challenge 
of  the  MoM  solution  is  determining  the  appropriate  weight,  or  coeffi¬ 
cient,  for  each  cell's  basis  function.  The  basis  function  coefficients 
are  found  by  replacing  the  original  boundary  condition  with  a  set  of 
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approximate  ones,  and  then  testing  each  basis  function  for  satisfaction 
of  the  approximated  boundary  condition  using  a  set  of  weights  or  testing 
functions  [21:221]. 

The  accuracy  and  computational  feasibility  of  the  MoM  solution  is 
determined  by  selection  of  the  appropriate  number  of  cells,  basis 
functions,  and  testing  functions.  Typically,  the  accuracy  of  the  MoM 
solution  increases  in  direct  proportion  to  the  number  of  the  cells  on 
the  body.  However,  the  number  of  cells  dictates  the  dimensions  of  the 
matrices  in  the  MoM  matrix  equation,  and  there  is  a  point  where  the 
increase  in  accuracy  is  offset  by  a  penalty  paid  in  increased  computer 
processing  time  and  memory  requirements.  Selecting  the  number  of  cells 
is  then  a  trade  between  acceptable  accuracy  and  demands  on  computer 
resources . 

The  accuracy  of  the  MoM  solution  is  also  affected  by  how  well  the 
basis  functions  satisfy  the  surface  current  boundary  conditions.  While 
complex  basis  functions  tend  to  improve  accuracy,  a  penalty  of  increased 
computation  difficulty  must  be  paid  [18:317], 

Finally,  selecting  the  appropriate  test  functions  can  also  affect 
the  MoM  solution  accuracy  and  complexity.  Since  the  test  functions  are 
integrated  against  the  basis  function  to  form  an  inner-product,  the 
complexity  of  the  test  function  also  dictates  the  complexity  of  the 
computations.  Secondly,  careful  selection  of  testing  functions  may  help 
to  avoid  integration  difficulties  if  the  fields  exhibit  singularities, 
such  as  those  of  the  fields  near  the  edges  of  a  conductor. 
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In  [18],  Mosig  addresses  the  issue  of  selecting  the  appropriate 
number  of  cells,  basis  functions,  and  testing  functions  for  a  radiating 
microstrip  patch.  In  this  research,  Mosig  used  several  schemes  to 
divide  an  arbitrarily  shaped  patch  while  measuring  the  computation 
requirements  and  accuracy  of  the  MoM  solution  for  different  combinations 
of  basis  and  testing  functions.  To  determine  the  best  choices,  Mosig 
considered  computation  speed  and  the  ability  to  predict  measured  values 
of  resonant  frequency,  quality  factors,  and  input  impedance  for  a  patch 
resonator  [18:322],  In  this  effort,  the  antenna  patch  was  divided  into 
rectangular  cells  referred  to  as  charge  cells,  with  two  adjacent  charge 
cells  forming  a  current  cell.  For  charge  cells,  Mosig  found  that  using 
cells  with  an  x-dimension  of  0.05X,  along  with  the  appropriate  basis  and 
testing  functions  provides  solution  accuracy  of  within  one  percent 
[18:320].  The  y-dimension  of  the  cells  does  not  affect  solution 
accuracy  very  much,  and  Mosig  reports  any  number  of  transverse  cells 
greater  than  3  should  prove  adequate  for  engineering  calculations 
[18:320].  The  basis  functions  of  choice  are  two-dimensional,  rooftop 
functions  defined  over  each  current  cell  [18:318].  The  x-directed 
surface  currents  and  y-directed  surface  currents  are  expanded  in 
independent  sets  of  rooftop  basis  functions.  With  the  rooftop  basis 
functions,  the  best  selection  for  testing  functions  are  a  unit  pulse 
defined  over  the  line  segment  connecting  the  centers  of  the  two  charge 
cells  comprising  a  single  current  cell  [18:318].  By  defining  the 
testing  functions  between  charge  cell  centers,  the  need  to  evaluate  the 
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fields  near  patch  edges  are  avoided  thus  the  field  singularities  at  the 
edges  will  not  adversely  effect  the  MoM  solution  [20:427], 


2.5  Summary  of  Current  Knowledge 

While  microstrip  antennas  offer  numerous  advantages  over  conven¬ 
tional  antennas,  their  applications  are  somewhat  restricted  due  to 
limited  bandwidth.  Pozar's  aperture  feed  technique  and  Sabban's  stacked 
patch  configuration  provide  methods  to  overcome  these  bandwidth  limita¬ 
tions  to  an  acceptable  degree.  Tsao  et  al  have  demonstrated  the  band¬ 
width  improvements  possible  by  combining  the  aperture  feed  and  stacked 
patch  configurations.  Nazar's  analysis  of  the  aperture  fed,  stacked 
patch  antenna  provides  the  Green's  functions  and  numerical  techniques 
necessary  to  develop  the  method  of  moments  solution  for  the  antenna 
currents.  Mosig's  research  in  applying  the  MoM  to  microstrip  radiators 
provides  the  means  to  maximize  the  accuracy  of  a  MoM  solution  while 
minimizing  computational  requirements.  Combining  the  results  from  Nazar 
and  Mosig  provides  a  proven  approach  for  employing  the  MoM  to  find  the 
electric  surface  currents  on  the  aperture  fed,  stacked  patch  microstrip 
antenna . 
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III.  Theory 


This  chapter  discusses  the  use  of  the  MPIE  Green's  functions  and 
the  method  of  moments  to  solve  for  the  surface  currents  of  the  aperture- 
fed,  stacked  patch  microstrip  antenna.  The  complete  antenna  model  is 
first  described,  and  the  equations  needed  to  find  the  unknown  surface 
currents  on  all  four  interfaces  of  the  antenna  are  explained.  The 
details  of  the  MoM  solution  are  described  and  the  equations  for  the  MoM 
matrix  elements  and  characteristics  of  the  resulting  matrix  equation  are 
presented. 

After  the  complete  model  of  the  antenna  is  derived,  some  of  the 
challenges  of  completing  the  solution  for  all  four  planes  of  the  antenna 
are  discussed.  Next,  the  full  model  is  modified  so  that  the  MoM 
solution  can  be  accomplished  independent  of  the  antenna  excitation 
technique  by  using  an  approximation  for  the  magnetic  fields  of  the 
antenna  aperture.  The  MoM  matrix  equation  for  the  resulting  3-plane 
solution  is  then  derived  from  the  results  for  the  full  antenna  model. 
Finally,  an  iterative  technique  for  solving  the  MoM  matrix  equation  and 
the  equations  needed  to  calculate  the  radiated  fields  of  the  antenna  are 
discussed. 
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Throughout  the  remainder  of  this  thesis,  vector  values  will  be 
designated  by  a  single  underline  and  a  dyadic  will  be  designated  using  a 
double  underline. 

3.1  The  Antenna  Model 

The  object  of  the  analysis  is  to  obtain  expressions  for  the 
electric  surface  currents  on  the  antenna.  To  begin  the  analysis,  the 
electric  and  magnetic  surface  currents  are  defined  as  shown  in  Figure 
3-1.  The  electric  surface  current  on  the  microstrip  feedline  produced 
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Figure  3-1  Equivalent  Antenna  Model  With  Incident  and  Induced  Currents, 
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by  the  applied  signal  is  defined  as  J.lnc ,  and  the  electric  surface 
current  on  the  feedline  resulting  from  electric  fields  internally 
scattered  by  the  antenna  is  defined  as  2f.  The  equivalence  principle  is 
used  to  replace  the  ground  plane  aperture  with  its  equivalent  magnetic 
surface  current,  Uj.  To  insure  continuity  of  the  tangential  electric 
field  through  the  aperture,  the  magnetic  surface  current  below  the 
ground  plane,  JJ“i.  is  defined  as  the  negative  of  that  above  the  ground 
plane,  ^lb1  [29:978].  The  electric  surface  current  vector  on  the  patch 
immediately  above  the  ground  plane  is  defined  as  J2  and  that  of  the 
uppermost  patch  is  J3.  The  space  below  the  ground  plane  is  designated 
region  a  (z  <  0)  and  that  above  the  ground  plane  is  region  b  (z  >  0) . 

The  antenna  coordinate  system's  origin  is  located  at  the  center  of  the 
aperture,  with  the  x-axis  extending  parallel  to  the  length  of  the 
microstrip  feedline. 

In  region  a,  the  feedline/dielectric  interface  is  numbered  2a  and 
the  dielectric/ground  plane  interface  is  numbered  la.  In  region  b,  the 
ground  plane/dielectric  interface  is  numbered  lb,  the  patch/dielectric 
interface  immediately  above  the  aperture  is  numbered  2b,  and  the 
uppermost  patch/dielectric  interface  is  numbered  3b.  While  only  the  x- 
components  of  the  surface  currents  are  shown  in  Figure  3-1,  this 
solution  considers  both  x-  and  y-directed  surface  currents. 

The  unknown  surface  currents  J.£,  J3,  and  are  found  by 

solving  the  four  MPIE  equations  that  describe  the  boundary  conditions 
for  the  tangential  magnetic  field  of  the  aperture  and  tangential 
electric  fields  on  the  feedline  and  patches  of  the  antenna.  These 
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equations  are  written  to  satisfy  the  boundary  conditions  of  Table  3-1 
and  are  solved  via  the  method  of  moments .  The  boundary  conditions  not 
shown  in  Table  3-1  are  incorporated  in  the  construction  of  the  Green's 
functions  for  the  antenna  [23:16]. 


Table  3-1  Antenna  Tangential  Field  Boundary  Conditions  [23] 


1) 

- 

0  on  the  Patch  of  Interface  2b 

2) 

- 

0  on  the  Patch  of  Interface  3b 

3) 

£t«n  _ 

0  on  the  feedline 

4) 

Ht,m  is 

continuous  through  the  aperture 

3.2  Tangential  Electric  Field  Equations 

Nazar's  analysis  provides  the  Green’s  function  necessary  to 
calculate  the  magnetic  vector  potential  and  electric  scalar  potential 
resulting  from  electric  surface  currents  and  electric  surface  charge 
density  via  a  superposition  integral  taken  over  the  source  surface.  For 
example,  to  find  the  magnetic  vector  and  electric  scalar  potentials  at 
interface  2b  resulting  from  the  electric  surface  currents  and  the 
electric  charge  density  on  interface  3b,  the  following  integral  expres¬ 
sions  are  used  [23:15] 

a!,u)  -J 

(3-1) 

V$23<£)  ‘  fS3bG$23(*)q3(£')ds' 
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where 

&b23(fi)  “  the  magnetic  vector  potential  produced  on  interface  2b 
by  surface  currents  on  interface  3b 

r  -  ip  -  p'  | 

S3b  »  the  surface  of  the  patch  on  interface  3b 

GbA23(R)  -  the  dyadic  Green's  function  used  to  calculate  the 

contribution  to  the  magnetic  vector  potential  on  interface  2b 

created  by  the  electric  surface  current  on  interface  3b 

Vbq23 (fi)  ”  the  electric  scalar  potential  on  interface  2b  created  by 

the  surface  charge  density  of  interface  3b 

Gbq23(R)  -  the  Green's  function  used  to  calculate  the  contribution 
to  the  electric  scalar  potential  at  interface  2b  created  by  the 
electric  charge  density  on  interface  3b 

q3C.fi')  “  the  surface  charge  density  on  the  patch  of  interface  3b 

The  dot  in  the  expression  for  the  magnetic  vector  potential 
indicates  the  vector  dot  product  of  the  dyadic  Green's  function  and  the 
electric  surface  current  vector.  The  first  numeral  in  each  subscript 
identifies  the  interface  where  the  field  is  observed  and  the  second 
numeral  identifies  the  interface  where  the  source  is  located.  The 
superscript  b  signifies  the  interfaces  of  interest  are  located  in 
region  b  of  Figure  3-1.  The  vectors  £  and  £'  are  the  position  vectors 
of  the  observer  and  source  points,  respectively.  Both  vectors  are  a 
function  of  the  cartesian  coordinate  pairs  (x.y)  and  (x',y')  or  the 
cylindrical  coordinate  pairs  ( p,<p )  and  respectively . 
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The  general  dyadic  Green's  function,  GAb(R)  is  defined 


Gb(£) 


fi£"(S)xx  G%xy(R)i  y 
Glyx(R)yx  Gl”(R)  yy 


(3-2) 


Using  the  time  dependence  of  exp (jut)  where  j  -  J-l,  the  vector 
and  scalar  potentials  are  then  used  to  calculate  the  tangential  electric 
fields : 


£  -  - ju>A  -  VV 
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The  tangential  electric  fields  generated  by  the  magnetic  surface 
currents  of  the  aperture  are  found  using  the  appropriate  Green's  func¬ 
tions  as  follows  [23:39]: 


era.)  -JSll  <*'>*' 


(3-6) 


where 

GbE2i(R,$’)  -  the  dyadic  Green's  functions  used  to  calculate  the 
contribution  to  the  electric  field  at  interface  2b  from  the 
magnetic  surface  current  of  interface  lb 
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The  variable  f  in  equation  (3-4)  is  the  angle  between  the  source 
direction  and  the  direction  to  observer  location,  and  is  defined  as 
follows : 


f 


r 

sin'1 

p  sin <p  -  p*  sin^  | 

L  *  J 

s  in'’  ^ 

p  cos <p  -  p'  COS(^  j 

L  *  J 

;  x-directed  sources 
;  y-directed  sources 
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On  interface  2b,  the  tangential  field  is  a  sum  of  contributions 
from  the  electric  surface  currents  on  both  patches  and  the  magnetic 
surface  current  of  the  aperture.  Using  equations  (3-1),  (3-3),  and 
(3-4),  the  integral  equation  for  the  tangential  electric  fields  at 
interface  2b  satisfying  the  first  boundary  condition  of  Table  3-1  is 
[23:40] 


-jwf  Gb(R)  «J2(£'  )ds'  -  Vb\  G\ZZ{R)  q2ie!  )  dS' 

Js2b“A22  Js2b 

-jwf  Cb  (R)  •  JA£!)dS'  -  Vtf  Gq23(J?)q3(£/)ds/  (3-6 

Js3b  — A23  “3  JS3b 

+  fc  *^1<£/) ds>  ~  0 

JSlb-E21 

vhere  the  superscript  t  on  the  del  operator  indicates  the  operation  is 
applied  only  to  the  transverse  observer  coordinates,  x  and  y  in  a 
cartesian  system  and  p  and  in  a  cylindrical  system. 
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Similarly,  the  integral  equation  for  the  tangential  electric 
fields  at  interface  3b  satisfying  the  second  boundary  condition  of 
Table  3-1  is  {23:40] 


->r  Gb(R)'  J2(£!)ds'  -  Vtf  G%32(R)q2(£')ds' 

Js2b  A32  Js2b 

->JS3b£A33(J?),^3(£/)dS'  ‘  7tJS3bGq33^)‘?3(£/)ds/  <3"7> 

*JSlb£L(K^>-a.^'>ds'  -° 


On  the  feedline  of  interface  2a,  the  tangential  electric  fields 
are  the  difference  between  the  tangential  fields  created  by  the  applied 
signal,  the  fields  created  by  the  reflected  current,  and  the  fields 
created  by  the  magnetic  currents  of  the  aperture.  The  equation  satisfy¬ 
ing  the  third  boundary  condition  of  Table  3-1  is  {23:41] 


-jujS2aG*22(R)-JiTlc(£!)ds'  -  VtJs2aG‘22(l?)9inc(£,)ds/ 

♦  jw  f  G‘(R)-JA£')ds'  *  Vlf  C‘2Z(R)  qt(£)ds'  (3-8) 

Js2a~AZZ  Js2a 

•  JSi 


The  final  boundary  condition  of  Table  3-1  requires  a  continuous 
tangential  magnetic  field  through  the  aperture.  Therefore,  the  tangen¬ 
tial  magnetic  fields  due  to  sources  above  the  ground  plane  must  equal 
those  fields  produced  by  sources  below  the  ground  plane.  To  calculate 
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the  magnetic  fields  produced  by  the  magnetic  surface  currents,  the 
appropriate  Green's  functions  are  first  integrated  against  the  magnetic 
surface  current  vector  and  magnetic  charge  density  to  obtain  the 
electric  vector  and  magnetic  scalar  potentials  [23:39] 


E(£)  -  Jg/  Gf(R)  •  MU')ds' 
V.(£)  -  J,,  <?.(*)  q.(£)  dsf 


(3-9) 


where 

F(p)  -  the  electric  vector  potential 

Gp(p,p')  ”  the  dyadic  Green's  function  used  to  calculate  the 
contribution  to  the  electric  vector  potential  from  the  magnetic 
surface  current 

Gm(p,p‘)  -  the  Green's  function  used  to  calculate  the  contribution 
to  the  magnetic  scalar  potential  from  the  magnetic  charge  density 
q m(p')  -  the  equivalent  magnetic  surface  charge  density 

Using  the  electric  vector  and  magnetic  scalar  potentials,  the 
magnetic  fields  are  then  calculated  from  the  relationship 


H  -  -JvF  -  VVa 


(3-10) 
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The  contribution  to  the  magnetic  fields  from  the  electric  surface 
current  vector  are  found  by  integrating  the  surface  current  vector 
against  the  appropriate  Green's  function  in  the  same  manner  as  equation 
(3-4). 

In  region  a,  the  magnetic  fields  are  created  by  the  electric 
surface  currents,  J.lnc  and  2f,  and  the  lover  magnetic  surface  current  of 
the  aperture,  Ui*-  In  region  b,  the  magnetic  fields  are  produced  by  the 
electric  surface  currents  on  both  patches,  Jj  and  Jj,  and  the  aperture's 
upper  magnetic  surface  current  of  the  aperture,  Mjb.  The  integral  equa¬ 
tion  satisfying  the  fourth  boundary  condition  of  Table  3-1  is  [23:41] 


G*  (R,S)'J.  (p')ds'  *  f  Ga  (£,0‘  JAe!)ds' 

2a=H12  5  '  JS2a‘'Hl2 

-  >JSiaG*u(l?)  .  h1<£/ )  ds'  -  Vt  JSiaG*u(i?)  qml(/  )  ds' 

-  jwJSib£^uW  *«!(£')  ds'  ♦  VtJs^G^n(J?)qml(£')ds' 


(3-11) 


Equations  (3-6),  (3-7),  (3-8),  and  (3-11)  are  the  four  equations 
that  are  solved  simultaneously  to  obtain  the  unknown  surface  currents  - 
J  f ,  ^ ,  J  ^ ,  and  • 
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3.3  Details  of  the  Method  of  Moments  Solution 


The  first  step  of  the  MoM  solution  is  the  division  of  the  patches, 
feedline,  and  aperture  into  rectangular  areas  referred  to  as  charge 
cells.  All  charge  cells  are  taken  to  be  the  same  size  and  have  dimen¬ 
sions  a  x  b,  see  Figure  3-2.  Research  completed  by  Mosig  [18:320]  has 
shown  the  most  suitable  charge  cell  dimensions  are  setting  a  -  0.05A  and 
choosing  b  such  that  there  are  3  or  more  charge  cells  across  the  plane 
in  the  y-direction.  Having  all  cells  the  same  size  is  not  an  a'  'olute 
requirement,  but  doing  so  reduces  the  complexity  and  computation  time 
for  the  solution  [20:425], 


Figure  3-2  Geometry  of  the  Antenna  Charge  and  Current  Cells 
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Two  adjacent  x-directed  (y-directed)  charge  cells  form  an  x- 
directed  (y-directed)  current  cell.  The  x-directed  line  segment 
connecting  the  center  of  the  jth  charge  cell  and  the  adjacent  charge 
cell  to  the  right  forms  an  x-directed  test  segment,  cxj  •  Similarly,  the 
y-directed  line  segment  connecting  the  center  of  the  kth  charge  cell 
and  the  adjacent  charge  cell  to  the  top  forms  a  y-directed  test  segment, 
Cyk.  This  arrangement  provides  an  overlap  of  current  cells  and  means 
the  number  of  current  cells  on  the  plane  depends  upon  the  number  of 
charge  cells.  The  exact  relationship  between  the  number  of  charge  cells 
and  current  cells  is  not  easily  determined  for  an  arbitrarily  shaped 
source.  However,  for  a  rectangular  patch  with  m  x  n  charge  cells,  the 
number  of  x-directed  current  cells  is  given  by  Nx  -  (m  -  l)n  and  the 
number  of  y-directed  cells  is  Ny  -  m(n  -  1)  [20:426],  Any  charge  cell 
not  on  the  perimeter  of  the  source  will  support  up  to  two  x-directed 
current  cells  and  two  y-directed  current  cells. 

The  center  of  test  segment  Cxj  and  the  jth  x-directed  current  cell 
has  position  vector  .  The  left  and  right  ends  of  CXj  have  position 
vectors  £"Xj  and  £+Xj ,  respectively.  The  position  vectors  for  the  test 
segment's  ends  are  related  to  the  center  position  vector  through  the 
relationships  [20:426] 


where  x  is  the  x-directed  unit  vector. 


(3-12) 
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Similarly,  the  y-directed  test  segment  Cy k  has  position  vectors 
£yk<  £~yk>  and  £%k-  The  vector  £~yk  identifies  the  position  of  the  center 
of  the  lower  charge  cell  in  the  current  cell  and  £*yk  identifies  the 
position  of  the  center  in  the  upper  cell.  The  relationship  between  the 
position  vectors  of  Cyk  is  given  by  equation  (3-12)  by  exchanging  b  **  a 
and  y  «*  x. 

3.3.1  Basis  and  Test  Functions .  To  obtain  the  most  general 
solution,  no  a  priori  assumptions  are  made  concerning  the  distribution 
of  currents  on  the  conductors  and  aperture.  This  stipulation  eliminates 
the  possibility  of  using  entire  domain  basis  functions  in  the  method  of 
moments  solution  [20:425].  Previous  comparisons  of  various  subdomain 
basis  functions  have  shown  that  rooftop  functions,  see  Figure  3-3,  can 
provide  an  acceptable  representation  for  the  surface  currents  while 
maintaining  relative  simplicity  in  the  MoM  impedance  element  calcula¬ 
tions  [ 18 : 320] . 

Each  x-directed  current  cell  supports  a  single  rooftop  function, 
Tx(p'),  described  by 


^(x7) 


1-i^LL;  |x'|<a,  \y'\<* 
0  ;  Otherwise 


(3-13) 


where 


x7  -  *•(*'-*'*) 

/  -  y-U' 


(3-14) 
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Figure  3-3  Surface  Current  Rooftop  Basis  Functions  and  Surface  Charge 
Pulse  Doublet  Basis  Functions  for  x-  and  y-directed  Current  Cells 

Similarly,  each  y-directed  current  cell  supports  a  rooftop  basis 
function  Ty(y').  The  expression  for  Ty(y')  is  obtained  by  interchanging 
a  **  b  and  x  «  y  in  equation  (3-13). 

Using  the  rooftop  basis  functions,  the  electric  currents  are 
expanded  in  the  finite  series 


J(£')  -  x  ^ 


T 

j-i 


t  ®4rT(/) 


(3-15) 


where 

Qj  -  the  unknown  surface  current  expansion  coefficients 
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All  four  unknown  surface  currents  are  expanded  in  the  same  manner. 
The  expansion  coefficients,  Qj ,  for  surface  current  sit  are  identified  by 
superscript  f,  the  coefficients  for  magnetic  surface  current  by 
superscript  1,  the  coefficients  of  £2  by  superscript  b2 ,  the  coeffi¬ 
cients  of  Jj  by  superscript  b3 ,  and  the  coefficients  Of  line  by  super¬ 
script  inc. 

As  discussed  by  Mosig  and  Gardiol  [20:426],  the  factors  of  1/a  and 
1/b  in  the  expressions  for  the  surface  currents  give  the  coefficients 
the  units  of  a  current.  In  fact,  the  coefficients  quantify  the  total 
current  flowing  across  the  adjacent  boundary  of  the  two  charge  cells 
that  form  each  current  cell  [20:426]. 

To  expand  the  electric  surface  charge  density,  the  surface  current 
expansion  is  substituted  into  the  continuity  equation 


V  •  J  +  Jug  -  0 


(3-16) 


and  the  surface  charge  density  expansion  is  [20:426] 


<?(£') 


Juab 


£  Qjn*<x')  +  £  Ojny(y': 

j-l 
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The  functions  n,(x')  -  -V1  •  Tx(x')  and  Ily(y')  -  -V*-  •  Ty(y')  are 
x-  and  y-directed  pulse  doublets  described  as  follows: 
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n x(p')  - 


n  y(p') 


1;  0 <x'<a,  | y<  \'<* 

-1;  -a<x'  <0,  |/  |<-| 
0;  Otherwise 

1;  |x'|<-|,  -b<yJ  <0 

-1;  |x'  |<*.  OiS/cb 

0;  Otherwise 

k 
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where  x'  and  y'  are  defined  by  equation  (3-14). 

The  charge  density  of  a  single  charge  cell  can  be  calculated  from 
the  expansion  coefficients  of  all  the  current  cells  sharing  the  charge 
cell.  For  example,  for  the  jth  charge  cell  identified  by  the  position 
vector  ,  see  Figure  3-2,  the  charge  is  calculated  using  [20:426] 


j  +  Qy.k*i 


Q 


y.k 


1 


(3-19) 


where  the  subscript  of  each  coefficient  is  taken  from  the  test  segments 
adjoining  the  adjacent  charge  cell  centers.  The  subscripts  j  and  j+1 
indicate  the  expansion  coefficients  belonging  to  the  current  cells 
formed  by  the  charge  cell  of  interest  and  the  charge  cell  to  the 
immediate  left,  subscript  j,  and  the  charge  cell  to  the  immediate  right, 
subscript  j+1 .  Similarly,  the  coefficient  k  identifies  the  coefficient 
for  the  current  cell  formed  by  the  charge  cell  of  interest  and  the 
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charge  cell  immediately  below,  and  the  subscript  k+1  identifies  the 
coefficient  for  the  current  cell  formed  by  the  charge  cell  of  interest 
and  the  charge  cell  immediately  above. 

Figure  3-3  shows  the  relationship  between  the  charge  cells, 
current  cells,  rooftop  functions,  and  pulse  doublets  for  both  x-  and  y- 
directed  current  cells. 

Using  the  same  rooftop  function  and  pulse  doublet  functions,  the 
expansions  for  the  magnetic  surface  currents  and  charge  densities  are 
obtained . 

With  rooftop  basis  functions,  a  practical  choice  for  test  func¬ 
tions  are  unidimensional,  rectangular  pulses  [ 18 : 318 ; 20 : 426 ] .  These 
pulse  test  functions  provide  a  simple  means  to  complete  the  required 
inner-products  for  the  MoM  solution  without  sacrificing  accuracy  of  the 
results.  The  test  functions  are  defined  only  along  the  paths  from  £~XJ 
to  £+XJ  on  x-directed  current  cells  and  £~ylc  to  £+yk  on  y-directed  current 
cells.  By  defining  the, test  segments  along  these  paths,  the  need  to 
compute  field  values  near  the  plane's  edges  is  eliminated,  thereby 
avoiding  field  singularities  that  adversely  effect  the  accuracy  of  the 
MoM  solution  [20:427], 

Now,  define  an  inner-product  for  the  x-directed  test  function, 
Wj(x)  and  any  general  vector  function  £(r)  such  that 


(3-20) 
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Similarly,  <wk(y),  £(p,p')>  is  defined  with  integration  along  the 
path  Cyk.  Since  the  test  functions  selected  for  this  solution  are  a 
unidimensional,  rectangular  pulse  of  amplitude  one,  the  inner  products 
reduce  to  integrations  of  the  function  Q( r)  along  the  appropriate  path. 

3.3.2  Converting  Integral  Equations  to  Matrix  Equations .  Follow¬ 
ing  the  standard  procedure  of  the  method  of  moments,  the  surface  current 
and  charge  density  expansions,  equations  (3-15)  and  (3-17),  are  substi¬ 
tuted  into  the  integral  equations  (3-6),  (3-7),  (3-8)  and  (3-11).  Next, 
an  inner-product  of  the  integral  equations  with  test  functions  are 
taken,  and  these  inner-products  form  the  elements  of  the  matrices  in  the 
MoM  matrix  equation. 

However,  before  proceeding  with  the  substitution  and  statement  of 
the  matrix  forms  of  the  integral  equations,  a  simplification  for  the 
field  contributions  from  the  surface  charges  should  be  discussed. 
Consider  the  fourth  term  of  equation  (3-6)  which  is  the  expression  for 
the  tangential  electric  field  on  interface  2b  created  by  the  surface 
charge  density  on  interface  2b: 


E\z{q2)  -  -V^Js^^22(^)<72(£/)ds/ 


(3-21) 


Note  that  this  expression  contains  a  gradient  operation  with  respect  to 
this  trsr.svsrss  obssri.,6r  cooirdi.p.a.t^s  x  and  y. 

Now  take  the  inner-product,  as  defined  by  equation  (3-20),  along 
the  path  Cxj  with  Wj(x)  -  1  and  £(r)  -  E1^^) 
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(3-22) 


(EjOO  .lL<9a>)  “ 


-L^*Vtjs3b^2W^^)dS/ 
Js3b  [fcxi 


q2(£/)ds/ 


The  fact  that  the  gradient  is  taken  with  respect  to  the  observer 
or  unprimed  coordinates  allows  an  exchange  in  the  order  of  operations  in 
the  inner-product.  Since  the  contour  Cjj  is  the  path  from  £-xj  to  £+xj , 
the  bracketed  integral  in  the  center  of  equation  (3-22)  can  be  simpli¬ 
fied  to  the  difference  of  the  Gbq22  evaluated  at  the  endpoints  of  CXj , 

£+xj  and  £"Xj. 

With  these  simplifications,  the  inner-product  of  the  fields 
created  by  the  surface  charge  and  the  test  functions  are 


(v,(x),£b2(q2))  -  -JsJ^22(£;.l£')-  ^22(£;.,£')]q2(/)ds'  (3-23) 


Using  the  inner  products  as  expressed  in  equation  ^3-23)  and  the 
expansion  for  the  surface  currents  and  charge  density  given  by  equations 
(3-15)  and  (3-17)  equation  (3-6)  describing  the  fields  at  interface  2b 
can  be  written 
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n2  n2 

E  E 

i-1  j-1 


h2  n3 

E  E 

i-l  j-l 


j;  -  ck<*o>  ]n<^>ds/ 


.  b2 


(3-24) 


b3 


iJ1{€k-^,£-(R-n-r(ij)<,s'H  ■ 0 


where 


Csi  “  kci  I 

Isis  N,x;  s-  1,  2,  3 

Csi  "  ^yi  I 

N„  <  i  <  Ns;  s  -  1.  2,  3 

N1(  N2,  N3 

-  the  total  number  of  current  cells 

on  interface 

interface 

lb,  interface  2b,  and  interface  3b, 

respectively 

tsj  “  as  ; 

1  <  j  <  Nsx,  s-1,  2,  3 

tSj  -  bs  ; 

Nsx  <  j  <  Ns,  s-1,  2,  3 

dl  -  x  dx 

;  Isis  Nsx,  s  -  1,  2,  3 

dl  -  v  dv 

;  Nsx  <  i  <  Ns ,  s  -  1,  2,  3 

j  -  index  of  the  current  cell  on  the  source  plane 
i  -  index  of  the  current  cell  on  the  observer  plane 
Rij  “  I  Pi  ”  P'  j! 

R+ij  “  |P+i  ~  P'  j  |  ,  R“ij  is  defined  similarly 
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£+xl  and  £~xi  -  the  position  vectors  of  the  endpoints  of  test 
segment  Cxi  as  defined  in  equation  (3-12) 

and 


is  ' 


I(T'j)  - 


n(7p 


*•(£'  -£'xj),0*J*N.x 
y  •(£'  -£'yJ).X.x<J*V. 


iTx(7j)  ;0<J*Ntx 

yTy(7j)  ;NX<J*N, 

nx(7j)  ;0<j<Ntx 

ny(7j)  ;Nsx<j<N, 


1.  2,  3 
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Equation  (3-24)  can  be  written  using  matrices  as  follows: 


(3-26) 


where  Z  indicates  impedance  matrices  ,'ith  units  of  ohms  per  unit  length, 
C  indicates  coupling  matrices  with  units  of  length-1  [23:45],  and  the 
superscripts  are  assigned  in  accordance  with  previous  conventions. 

The  same  substitutions  and  simplifications  for  equation  (3-6), 
which  describes  the  fields  on  interface  3b,  produces 
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n3  n2 

E  E 

i- 1  j- 1 


A!±  f  dl  •  f  Gb(Ri.)  .  r(7,)ds/ 
^2J  JC3i -  Js2j“A32  -WJ/ 


b2 


»3  »3 

E  E 

i-i  j-i 


K2  N! 


Jj±  f  dl-  f  Gb  (J?n)  •  Tii^ds' 
t3jJc3i—  Js3j-A33V  -WJ' 


(3-27) 


r«$3 


The  equivalent  matrix  equivalent  for  equation  (3-27)  is 


>*»][»«].  [zf»][.“]-[cS‘][.5]-o  <3-28> 


Before  substituting  the  expansion  for  the  surface  current  and 
charge  density  into  the  equation  for  the  fields  of  interface  2a,  the 
terms  are  rearranged  so  that  the  known  quantity  Jinc  is  on  the  right  of 
the  equal  sign.  Using  the  subscript  f  to  identify  variables  associated 
with  the  microstrip  feedline,  equation  (3-8)  can  be  written 


3-22 


Nf  N  f 

E  E 

i-i  j-i 


Jj±  f  dl-  f  Gb  (Ji..)  •  r(7,)ds/ 

tfjJCfi—  Js£j-A22V  ~Wj' 

JSf  Jc^CKiV  -  C$22(*i'j)  ]n(7j)  ds' 


jwaf 


(3-29) 


E  E 

i-X  j-1 


_I_  f  di •  f  g*  (jji j  •  r(7,)ds/ 

tfjJc£1—  Js£j-A22V  -KJy 


inc 


which  produces  the  matrix  equation 
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Finally,  equation  (3-11)  describing  the  continuity  of  the  aperture 
magnetic  fields  becomes 
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These  four  matrix  equations  are  solved  simultaneously  to  form  the 
MoM  matrix  equation.  The  impedance,  coupling,  and  admittance  matrices 
of  equations  (3-26),  (3-28),  (3-30),  and  (3-32)  are  sub-matrices  in  the 
MoM  solution  matrix.  The  complete  MoM  matrix  equation  is 


to33]  to33] 

-to31] 

to] 

to) 

[0] 

fer] 

[afr] 

-to31] 

[0] 

to] 

[0] 

[,,bl2]  tobl3 

ru  J  . 

to] 

-to13] 

l-i] 

to13] 

[o] 

[o] 

to!,l 

[*?f  ]  . 

k] 

to33] 

(3-34) 


3.3.3  Properties  of  the  MoM  Solution  Matrix.  The  diagonal  sub¬ 
matrices  in  the  MoM  solution  matrix  of  equation  (3-34)  characterize  the 
contributions  to  the  tangential  fields  on  an  interface  created  by 
surface  currents  located  on  the  same  interface.  The  diagonal  terms  of 
each  of  these  sub-matrices  characterize  the  fields  of  a  current  cell 
caused  by  the  currents  in  that  same  cell.  These  self-impedance  are  the 
strongest  and  thus,  each  sub-matrix  and  the  solution  matrix  is  diagonal¬ 
ly  dominant  [23:47].  While  the  MoM  impedance  matrix  is  considerably 
large,  it  can  be  inverted  using  standard  Gaussian  elimination  [20:427]. 
However,  a  solution  using  an  iterative  scheme  offers  a  savings  in 
computer  processing  time  [17:271]. 

The  position  of  each  sub-matrix  of  the  solution  matrix  determines 
its  role  in  describing  the  tangential  fields.  The  row  index  identifies 
the  observer  location  with  rows  1  through  4  corresponding  to  interfaces 
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2b,  3b,  la  and  b,  and  2a,  respectively.  The  column  index  identifies  the 
source  location  with  the  same  order  as  the  row  index.  Thus,  the  sub¬ 
matrix  located  in  row  1  and  column  2  characterizes  the  fields  observed 
on  interface  2b  resulting  from  sources  located  on  interface  3b.  The 
zero  sub-matrices  identify  locations  in  the  antenna  where  sources  do  not 
contribute  to  observed  fields.  For  instance,  the  zero  matrix  in  row  2, 
column  4  signifies  the  sources  of  interface  2a  do  not  contribute  to  the 
fields  observed  at  interface  3b. 

The  number  of  rows  in  each  sub-matrix  is  determined  by  the  total 
number  of  current  cells  on  the  observer  plane,  and  the  number  of  columns 
is  determined  by  the  total  number  of  current  cells  in  the  source  plane. 
Each  diagonal  sub-matrix  will  be  a  square  matrix  because  the  observer 
and  source  planes  are  the  same  for  these  matrices.  Since  each  source 
and  observer  plane  is  represented  in  both  dimensions  of  the  solution 
matrix,  the  solution  matrix  is  a  square  matrix  with  its  dimension  in 
either  direction  being  equal  to  the  sum  of  the  current  cells  on  all  four 
interfaces  of  the  antenna. 


3.4  Challenges  of  the  MoM  Solution 

Completing  the  MoM  solution  for  the  entire  antenna  is  an  arduous 
process.  Before  any  calculations  on  the  MoM  solution  are  possible,  the 
Sommerfeld  integrals  defining  the  Green's  functions  must  be  evaluated, 
and  the  expansion  coefficients  for  jiine  must  be  calculated. 
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3.4.2  Evaluating  the  Antenna's  Green's  Functions .  Due  to  the 
numerical  complexity  of  evaluating  the  Sommerfeld  integrals,  the  Green's 
functions  for  the  antenna  will  not  be  explicitly  evaluated  in  the 
software  written  to  complete  the  MoM  solution. 

As  shown  by  Nazar  [23:74-84)  and  discussed  by  Mosig  and  Gardiol 
[21:224-226],  the  Green's  functions  are  dependent  only  on  kR,  where  k  is 
the  wavenumber  of  the  input  signal.  The  Green's  function  can  therefore 
be  evaluated  at  a  finite  number  of  points  within  a  range  of  kR  and 
interpolation  can  be  used  to  find  the  Green's  .functions  anywhere  within 
the  interval  [ 19 : 429  ;  19  :  224]  . 

The  approach  used  here  will  be  to  employ  polynomial  representa¬ 
tions  of  the  general  Green's  functions  as  presented  by  Erwert  [9].  When 
the  asymptotic  form  of  a  Green's  function  for  R  -*  0  is  needed  a  combina¬ 
tion  of  interpolation  from  tabulated  data  and  an  analytic  term  will  be 
used.  The  asymptotic  approximations  are  discussed  in  Sections  3. 4. 1.1 
and  3. 4. 1.2.  The  approximations  of  the  Green’s  functions  are  either 
sixth-order  polynomials  or  a  combination  of  an  eighth  or  tenth-order 
polynomials.  Because  the  Green's  functions  are  also  dependent  on  the 
applied  signal's  wavenumber,  different  polynomials  must  be  used  as  the 
input  frequency  is  varied. 

3. 4. 2. 1  Asymptotic  Green's  Functions  for  R  -*  0 .  For  the 
case  when  source  and  observer  locations  are  coplanar  and  on  the  same 
current  cell,  the  value  of  R  will  approach  zero.  Since  the  Sommerfeld 
integrals  defining  the  Green's  functions  contain  Hankel  functions  of  the 
second  kind,  the  integral  will  be  undefined  when  R  -*  0.  To  overcome  the 
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singularity  at  zero,  Nazar  derived  an  asymptotic  form  for  the  Sommerfeld 
integrals  of  interest  -  Gb1*^,  GbA33,  Gbq22  and  Gbq33.  The  asymptotic 
expression  involves  two  integrals  of  Bessel  functions  over  finite 
intervals  and  a  constant  term  proportional  to  R-1.  For  example,  the 
asymptotic  form  of  the  Green's  functions  describing  x-directed  fields  on 
interface  2b  created  by  x-directed  surface  currents  on  interface  2b  is 
given  by  [23:59] 


/-.bxx 

°A22 


(R) 


^ib 

~Tk 


r *c  Jo(  ^R)  ^ 

u2bcosh(u2b  (b2b-blh)) 

J°  Db(A) 

*  Mb23u3bs^n^(u2b  (^2b-^lb  )  ) 

1  *t*b23  Jp(x)  dx  +  _ l*Mb23 

1+Mbl3  +  ^bl2+Mb23  J°  R  R(l*  /^bl3+/ibl2+Mb23) 


(3-35) 


The  variables  in  the  first  and  third  terms  of  equation  (3-35)  are 
defined  on  the  first  page  of  Appendix  A. 

The  integrals  of  equation  (3-35)  are  now  over  finite  intervals, 
are  therefore  finite  for  all  R,  and  can  be  evaluated  numerically.  The 
thiru  term  of  the  approximation  must  be  integrated  analytically  over  the 
surface  of  the  source  for  use  in  the  method  of  moments  solution.  The 
evaluation  of  the  integral  is  the  subj  t  of  the  following  section. 

3. 4. 1. 2  Evaluating  the  Analytic  Term  in  the  Asymptotic 
Green's  Function  Expression.  With  the  rooftop  basis  functions,  the 
integral  defining  the  surface  integration  over  the  jth  x-directed 
current  cell  for  the  R_1  term  in  the  asymptotic  expression  is  given  by 
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(3-36) 


HR) 


ds '  dx 


where 

Xj,  yj  -  the  x-  and  y-components  of  the  cell's  position  vector, 
ds'  -  dy'  dx' 

R  -  [(x  -  x')2  +  (yj  -  y')2]* 


Completing  the  surface  integration  produces  the  expression 

r 

Q  2  2 

csch*1  (3r)  +  c1sinh'1(3r)  -  ^  csch'1  (r)  +  c2sinh~1  (r) 

2” 

+  ^(3r2*lp-  l(r2  +  l)'*-|(3r2  +  lp  +  J£(r2  +  lp 


where 


54a2  -  7b2  +  18ab 
TZTTE 

bz -  42a 2  -  1 8 a  b 
24  ab 


(3-38) 


The  result  of  the  surface  integration  over  a  y-directed  current 
cell  is  obtained  by  interchanging  a  «  b  in  equations  (3-37)  and  (3-38). 
The  details  of  completing  the  integration  of  equation  (3-36)  are  given 
in  Appendix  B. 


3-29 


Using  the  result  of  equation  (3-37),  the  asymptotic  value  of  the 
Green's  function  can  be  obtained  as  follows: 

1)  Using  tabulated  data  and  linear  interpolation,  find  the 
contribution  to  the  Green's  function  from  the  integral  of 
equation  (3-35). 

2)  Using  numerical  integration,  complete  the  surface  inte¬ 
grations  over  the  source  cell  for  the  product  of  the  inte¬ 
grals  of  equation  (3-35)  and  the  rooftop  basis  function. 

3)  Add  the  contribution  of  equation  (3-38),  multiplied  by 
the  appropriate  constant,  to  the  result  of  step  #2.  The 
result  is  the  integration  of  the  product  of  the  Green’s 
function  and  the  rooftop  basis  function  over  the  source 
cell's  surface. 

3.4.2  Modeling  the  Microstrip  Feedline.  To  complete  the  MoM 
solution  as  described  in  section  3.3,  the  expansion  coefficients  for  the 
surface  current  Jinc  must  be  known.  Finding  the  coefficients  requires 
completion  of  another  MoM  solution  using  the  same  basis  and  test 
functions  as  those  of  the  antenna's  solution.  The  exact  excitation 
model  would  depend  upon  the  means  used  to  excite  the  microstrip  feed¬ 
line.  One  common  excitation  is  the  use  of  a  coaxial  probe,  and  several 
means  to  model  this  excitation  have  been  described  by  various  authors 
[ 5 : 613 ; 19 :431-436 ; 20 :428 ; 21 : 229-231 ]  with  one  especially  detailed  model 
developed  by  Hall  and  Mosig  for  microstrip  antennas  with  electrically 
thick  substrates  [10:370-375],  For  a  general  mircostrip  feedline, 
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Kobayashi  presents  closed-form  expressions  for  the  normalized  longitudi¬ 
nal  and  transverse  feedline  current,  distributions  [20], 

Since  the  excitation  model  can  cake  on  such  wide  variations,  and 
the  purpose  of  this  research  is  to  obtain  a  more  general  solution  for 
the  antenna,  the  MoM  solution  of  section  3.3  will  be  changed  to  elimi¬ 
nate  dependence  upon  the  antenna  excitation  technique.  The  modified 
version  of  the  MoM  solution  is  explained  in  the  following  section. 


3.5  The  3 -Plane  Model  -  A  SoluCion  Independent  of  the  Excitation  Source 

To  prevent  the  MoM  solution  from  being  dependent  upon  the  excita¬ 
tion  model,  the  antenna  model  can  be  generalized  by  assuming  a  distri¬ 
bution  for  the  magnetic  currents  on  the  aperture.  By  making  this 
assumption,  the  feedline's  role  in  the  solution  is  eliminated,  and  the 
solution  now  depends  only  upon  the  interactions  of  the  aperture  and 
patches . 

3.5.1  Approximations  for  the  Aperture  Fields.  As  discussed  by 
Pozar  [26:50],  the  best  coupling  of  energy  from  the  aperture  to  the 
driven  patch  is  obtained  when  the  aperture  is  a  long,  thin  rectangular 
slot.  With  the  geometry  of  this  analysis,  the  long,  thin  rectangular 
slot  results  in  the  y-dimension  of  the  aperture  being  considerably 
larger  than  the  x-dimension.  Given  this  configuration,  the  electric 
fields  for  the  aperture  will  be  dominated  by  an  x-directed  component 
resulting  in  a  y-directed  equivalent  magnetic  current  on  the  aperture. 
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Since  the  aperture  is  electrically  small,  it  is  reasonable  to  assume  the 
electric  field  of  the  aperture  can  be  represented  by  a  single  piecewise 
sinusoidal  function  (27 : 1443 ; 27 : 978 ] 


(3-39) 


where 

W  -  the  x-dimension  of  the  aperture 
L  -  the  y-dimension  of  the  aperture 
ke  -  the  effective  wavenumber  of  the  piecewise  sine 

The  effective  wavenumber  is  best  represented  as  the  average  of  the 
wavenumbers  in  the  regions  adjacent  to  the  aperture  [ 26 : 1440  ;  29  : 979 ]  . 
Thus,  the  effective  wavenumber  of  the  piecewise  mode  is  given  by 


k,  - 


£la*€lb 
c  r  e  r 


(3-40) 


where 

k0  -  free  space  wavenumber 

€rla  -  the  relative  permittivity  of  the  dielectric  between  inter¬ 
faces  la  and  2a  of  the  antenna 
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€rlb  -  the  relative  permittivity  of  the  dielectric  between  inter¬ 


faces  lb  and  2b  of  the  antenna 


With  this  approximation  for  the  aperture  electric  fields,  the 
magnetic  currents  can  be  calculated  as  follows 


sin 


Wb  -  n  x  £b  «  z  x  xE\ x  -  y 


W  sin 


.fe-m 
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By  approximating  the  aperture  magnetic  currents  in  this  manner, 
the  solution  of  the  antenna  becomes  a  matter  of  satisfying  the  first  two 
boundary  conditions  in  Table  3-1. 

3.5.2  Modified  Matrix  Equation.  As  discussed  previously,  on  each 
patch  the  tangential  electric  fields  are  the  sum  of  contributions  from 
the  magnetic  current  of  the  aperture,  the  surface  currents  on  the  patch 
itself,  and  the  surface  currents  of  the  remaining  patch.  Since  the 
magnetic  currents  of  the  aperture  are  now  a  known  variable,  the  fields 
created  by  the  magnetic  currents  become  the  forcing  function  in  the 
field  equations.  Thus,  the  matrix  equation  satisfying  the  first  two 
boundary  conditions  of  Table  3-1  is 


3-33 


(3-42) 
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where  all  the  elements  on  the  left  of  the  equation  are  the  same  as  those 
previously  specified  in  equations  (3-27)  and  (3-31). 

The  inner-products  comprising  the  elements  of  the  column  matrix  to 
the  right  of  the  equals  sign  in  equation  (3-42)  are  given  by 


(3-43) 


Equation  (3-42)  is  the  method  of  moments  matrix  equation  that  will 
be  solved  in  this  thesis. 


3.6  Solving  the  MoM  Matrix  Equation 

The  solution  matrix  of  equation  (3-34)  and  impedance  matrix  of 
equation  (3-42)  both  have  large  dimensions.  For  example,  consider  an 
antenna  with  dimensions  shown  in  Table  2.  Using  the  charge  cell 
dimension  of  a  -  0.05A  and  choosing  five  y-directed  charge  cells  across 
each  plane  of  the  antenna,  the  dimensions  of  the  square  solution  matrix 
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of  equation  (3-34)  would  be  448,  and  with  the  3-plane  model,  the  square 
impedance  matrix  would  have  dimensions  of  170,  Solving  the  MoM  matrix 


Table  3-2  An  Example  Antenna  with  its  Associated  MoM  Segmentation 


Antenna  Element 

x-dimensions 

(wavelengths) 

x-directed 
current  cells 

y-directed 
current  cells 

Feedline 

1.5 

145 

120 

Aperture 

0.1 

5 

8 

Lower  Patch 

0.5 

45 

40 

Upper  Patch 

0.5 

45 

40 

equation  of  the  3-plane  model  using  Gaussian  elimination  to  invert  the 
impedance  matrix  is  feasible  [20:427].  However,  for  the  larger  matrix 
of  the  complete  antenna  model,  an  iterative  solution  scheme  might  prove 
to  be  more  advantageous . 

3.6.1  Solving  the  Matrix  Equation  via  Iteration .  As  discussed  by 
Moore  and  Pizar  [17:270],  iterative  schemes  offer  economies  of  computer 
storage  and  processing  time,  as  compared  to  Gaussian  elimination  or 
other  factorization  approaches,  when  the  matrices  in  the  MoM  solution 
are  very  large.  Since  the  square  impedance  matrices  of  the  three  or 
four  plane  MoM  solution  may  have  very  large  dimensions  for  certain 
geometries,  development  of  an  iterative  scheme  for  solving  the  MoM 
matrix  equation  is  worthwhile. 

3.6.2  Basic  Iterative  Method.  For  the  matrix  equation  ZI  -  V, 
the  matrix  Z  may  be  decomposed  such  that 
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Z  -  2d  +  2} 


(3-44) 


where  Z0  is  more  easily  inverted  than  Z.  By  substituting  equation 
(3-44)  into  the  original  matrix  equation,  the  equation  may  be  rewritten 


Z0I  -  V  -  za 


(3-45) 


Based  on  equation  (3-44),  an  iterative  process  can  be  defined  such 
that  [16:272] 


2o j (k ♦  l )  _  ]/  _  ;  k-0,1,2, .  (3-46) 

Moore  and  Pizer  discuss  several  iterative  schemes  based  on  a 
priori  knowledge  of  various  characteristics  of  the  matrix  system  being 
solved.  Earlier  discussions  of  the  MoM  matrix  equation  and  characteris¬ 
tics  of  the  Z  matrix  have  emphasized  the  fact  that  the  MoM  solution  and 
impedance  matrices  are  diagonally  dominant,  and  one  approach  highlighted 
by  Moore  and  Pizer  is  well  suited  to  a  system  with  this  characteristic. 

3.6.3  An  Iterative  Approach  Suited  to  the  MoM  Solution .  The 
iterative  scheme  best  suited  to  the  MoM  solution  of  the  complete  antenna 
model  was  proposed  by  Ferguson  et  al.  [9:230-235]  and  begins  with  the 
decomposition  of  Z  such  that 
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Z  -  I  +  B  +  U 


(3-47) 


where  B  is  a  banded  matrix  formed  by  elements  on  either  side  of  the 
diagonal  of  Z.  The  matrices  L  and  U  are  the  lower  and  upper  triangular 
matrices,  respectively,  left  when  the  elements  of  B  are  replaced  by  zero 
in  Z. 

The  key  to  this  scheme,  making  it  ideally  suited  to  the  solution 
of  the  aperture-fed,  stacked  patch  antenna,  is  that  the  matrix  B  must 
contain  the  largest  elements  of  Z  in  order  for  the  system  to  converge. 
Since  the  Z  matrix  of  the  MoM  solution  is  diagonally  dominant,  this 
condition  will  always  be  satisfied.  Also,  since  the  largest  elements  of 
the  MoM  matrix  are  already  on  the  diagonal,  no  row  pivoting  is  required 
during  decomposition  which  results  in  the  most  efficient  iterative 
solution  scheme  [9:230], 

Using  this  representation  for  Z,  the  iterative  scheme  becomes 


(3-48) 


The  iteration  is  repeated  until  a  previously  selected  stopping 
criterion  is  satisfied.  One  criterion  presented  by  Moore  and  Pizer 
[17:276]  is  based  on  the  relative  change  in  the  norms  of  successive 
solutions.  Denote  the  present  iterated  solution  with  a  superscript 
(k+1)  and  the  last  iteration  with  superscript  (k) ,  and  suppose  the 
vector  X  contains  N  elements.  The  first  step  in  determining  satisfac- 
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tion  of  the  stopping  criterion  is  to  form  the  difference  vector  6I(k+1) 
where 


5J( k.i)  _  J(k-1)  _  J(k)  (3-49) 

where  the  ich  element  of  the  vector  SXtk+1>  is  the  difference  between  the 
ith  element  of  I(k+1)  and  I(k).  Next,  the  norm  of  6X<k+1>  is  calculated 

|51<k*i>|2-  E  (3-50) 

i-1 

where  the  *  indicates  complex  conjugate.  Similarly,  the  norm  of  itk)  is 
found.  The  stopping  criterion  is  satisfied  when 

I  6I(k)  |  <  e  |  I(k*1)  |  (3-51) 


where  e  is  selected  based  on  the  required  accuracy  of  the  solution. 

Ferguson  et  al .  report  processor  time  saving  of  5  to  10  times,  and 
up  to  23  times,  relative  to  Gaussian  elimination.  However,  this 
particular  scheme  does  not  allow  any  computer  memory  savings  unless  some 
algorithm  for  disc  storage  of  the  matrices  in  the  scheme  is  employed. 
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3.7  Calculating  the  Fields  of  the  Antenna 


While  the  calculations  of  the  fields  on  the  antenna  interfaces 
will  not  be  a  part  of  this  effort,  the  software  developed  to  complete 
the  MoM  solution  can  be  used  to  do  so.  This  section  outlines  how  the 
fields  on  the  antenna  interface  and  the  far-field  patterns  can  be 
calculated  using  the  results  of  the  MoM  solution. 

3.7.1  Calculating  the  fields  on  the  Antenna  Interfaces .  After 
the  surface  current  and  charge  density  coefficients  are  known,  all  the 
sources  on  the  antenna  are  defined,  and  the  fields  at  any  location  on 
either  of  the  three  interfaces  of  the  3-plane  antenna  model  can  be 
calculated  using  many  of  the  same  Green's  function  da-a  as  used  in  the 
MoM  solution.  The  only  additional  data  needed  is  the  analytic  expres¬ 
sion  for  the  divergence  of  the  surface  charge  Green's  functions  so  that 
the  surface  charge  field  contributions  of  equations  (3-6)  and  (3-7)  can 
be  calculated. 

Once  the  divergence  of  the  surface  charge  Green's  functions  is 
known,  the  subroutines  of  the  MoM  software  can  be  used  to  complete  the 
integrations  needed  to  find  the  fields  at  any  point  on  the  antenna's 
interfaces . 

One  reason  to  calculate  the  fields  on  the  interfaces  is  that  doing 
so  would  provide  a  means  to  validate  the  MoM  solution.  If  the  solution 
is  correct,  the  fields  at  any  location  on  either  of  the  two  patches 
should  be  approximately  equal  to  zero.  Another  reason  to  calculate  the 
fields  would  be  to  calculate  the  radiation  patterns  of  the  antenna.  One 
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approach  to  finding  the  antenna  patterns  is  presented  in  the  following 
section . 

3.7.2  Calculating  the  Antenna  Far-Field  Patterns.  As  mentioned 
in  Section  1.1,  radiation  from  microstrip  antennas  is  produced  by  the 
fringing  fields  along  the  edges  of  the  microstrip  patch.  In  the 
aperture-fed,  stacked  patch  microstrip  antenna  the  fringing  fields  of 
the  patch  on  interface  3b  produce  the  radiation.  One  approach  to 
finding  the  antenna  field  patterns  is  to  model  the  fringing  fields  of 
the  patch  with  equivalent  radiating  apertures.  The  far-field  patterns 
of  the  antenna  can  then  be  calculated  using  theory  for  radiating 
apertures . 

3. 7. 2.1  An  Aperture  Model  for  the  Antenna.  Hammer  et  al .  have 
shown  the  validity  of  using  an  aperture  model  to  calculate  the  radiated 
fields  of  a  microstrip  antennas  [12].  The  first  step  of  the  modelling 
is  the  construction  of  four  apertures  along  the  edges  of  the  antenna's 
radiating  element.  Next,  the  tangential  fields  on  the  interface  of  the 
radiating  patch  are  calculated,  and  an  equivalent  electric  field 
distribution  is  found  such  that  the  following  relationship  is  satisfied 


J"  IIt | *dl  -  Ea* 2a  (3-52) 

where 

Et  -  the  tangential  electric  field  value 
Ea  -  the  equivalent  field  distribution  value 
2a  -  the  width  of  the  radiation  aperture 
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Actually,  the  value  of  £t  decreases  to  zero  at  a  finite  distance 
from  the  patch  edge,  and  the  upper  limit  of  the  integral  in  equation 
(3-52)  can  be  replaced  with  a  distance  L  where,  in  general,  L  is  much 
less  than  the  value  of  the  free-space  wavelength  (12:267). 

In  [12],  a  closed-resonator  model  is  used  to  approximate  the 
tangential  field  magnitude.  In  this  analysis,  the  Green's  functions  and 
the  electric  and  magnetic  surface  currents  permit  the  calculation  of  the 
tangential  fields.  Once  the  tangential  fields  have  been  calculated,  the 
actual  x-component  and  y-component  of  the  fields  can  be  replaced  with  an 
equivalent  distribution  such  that  the  relation  of  equation  (3-52)  is 
satisfied  for  the  x-directed  and  y-directed  fields.  This  procedure 
produces  an  aperture  model  for  the  antenna  as  illustrated  in  Figure  3-4. 
Since  the  radiating  apertures  of  the  antenna  model  have  constant, 
uniform  field  distributions ,  the  far-field  patterns  are  readily  calcu¬ 
lated. 

3. 7. 2. 2  Far-Field  Pattern  Equations.  The  far-field  pattern  is 
the  sum  of  the  pattern  produced  by  the  radiation  of  each  aperture.  To 
begin,  consider  aperture  1  of  Figure  3-4  with  equivalent  tangential 
fields  represented  by  the  following  expression: 


iti  -  +  Eyly 


(3-53) 
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Figure  3-4  Aperture  Model  for  Calculating  Far-Field  Radiation  Patterns 

The  equivalent  field  principle  can  be  used  to  obtain  an  expression 
for  the  magnetic  surface  currents  of  the  aperture  [2:460J: 


—ai  ~  2n  x  Etl  -  2 Eylx  -  2Exly 


(3-54) 


Since  the  magnetic  surface  current  can  be  expressed  in  terms  of 
the  components  of  the  tangential  fields,  the  components  of  the  far-field 
pattern  of  the  aperture  can  be  found  using  the  relationships  [2:455-458] 


E91  -  [ -Eyl  sin<p  -  Exl  cosip ]  eJkQr/  C05^  ds' 

E  ifi  "  [ £yi  cos  8  cos  i p  -  Exl  cos  8  sin<p  ]  Jg/  eJk£,r/  cos*  ds; 


(3-55) 
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where 


S'  -  the  aperture  surface 

i>  -  x'  cos 6  cos<p  +  y'  cos 6  simp 


and  the  constants  in  the  expression  have  been  neglected  since  only  the 
field  pattern  equations  are  desired. 
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IV.  Method  of  Moments  Software  Description 


This  chapter  presents  the  details  of  the  FORTRAN  software  devel¬ 
oped  to  complete  the  MoM  solution  for  the  3-plane  model  of  the  aperture- 
fed,  stacked  patch  microstrip  antenna.  The  chapter  begins  by  deriving 
the  numerical  integration  equations  needed  to  perform  the  various 
integrations  of  the  MoM  solution.  Next,  the  algorithms  and  logic  of  the 
subroutines  of  the  MoM  solution  software  are  described. 


4.1  Numerical  Integration  Equations 

The  major  task  of  completing  the  MoM  solution  for  the  antenna  is 
performing  the  integrations  needed  to  evaluate  the  elements  of  the  MoM 
matrix  equation.  In  the  interest  of  accuracy  and  economy,  Gaussian 
quadrature  is  used  to  complete  all  of  these  integrations. 

To  begin,  the  number  of  integration  nodes  to  be  used  in  the 
Gaussian  quadrature  must  be  determined.  The  selection  of  the  number  of 
nodes  is  based  on  the  numerical  behavior  of  the  integrand  within  the 
interval  of  integration.  If  the  integrand  is  oscillatory  within  the 
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interval,  a  larger  number  of  nodes  is  required  to  accurately  evaluate 
the  integral.  Therefore,  the  numerical  stability  of  the  integrands  of 
the  matrix  elements  must  be  investigated  before  a  selection  of  the 
number  of  nodes  is  made . 

Since  the  integrands  are  products  of  a  Green's  function  and  either 
a  linear  rooftop  function  or  the  constant  II  function,  the  numerical 
characteristics  of  the  Green's  function  will  dominate  the  integrand's 
behavior.  For  the  surface  and  line  integrals  in  the  MoM  solution,  the 
arguments  of  the  Green's  functions  are  rather  small.  Based  on  data 
presented  by  Nazar  [23:74-84],  under  these  conditions  the  real  and 
imaginary  parts  of  the  Green's  functions  are  well  behaved.  Therefore  a 
relatively  small  number  of  quadrature  nodes  should  suffice.  To  simplify 
software  development  and  use,  the  same  number  of  quadrature  nodes  will 
be  used  for  all  integrations .  For  the  remainder  of  this  discussion, 
five  node  quadrature  will  be  assumed. 

To  derive  the  numerical  integration  equations,  the  integrations 
for  the  impedance  sub-matrix  describing  the  fields  on  interface  2b 
produced  by  the  surface  currents  and  charge  density  of  interface  3b  will 
be  considered: 
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This  expression  is  typical  of  the  expressions  defining  the 
remaining  sub-matrices  of  the  MoM  impedance  matrix  and  the  numerical 
integration  equations  for  this  sub-matrix  can  easily  be  adapted  for 
those  sub-matrices. 

To  begin,  the  dot  product  operations  in  equations  (4-1)  must  be 
completed.  First,  the  dot  product  of  the  dyadic  Green's  function  and 
the  rooftop  basis  functions  in  the  source  surface  integration  will  be 
simplified.  Recall  for  1  <  j  s  N3x,  the  rooftop  vector  functions  are  x- 
directed,  and  for  N3x  <  j  <  N3,  the  rooftop  vect  r  functions  are  y- 
directed.  Therefore,  the  dot  product  between  the  dyadic  Green's 
function  and  the  rooftop  vector  functions  produces  the  following 
vectors : 


Z<73j) 


.[xGAb^(^J)+yGAb^(^J)]T7(73j)  !  l*j^3x 
[xG^fJCR^)  +  y Gb2y3y(Rij)]ry(73J)  ;  N3x<j<N3 


(4-2) 


Next,  consider  the  line  integration  along  the  observer  current 
cell's  test  segment  dotted  with  the  surface  integration  over  the  source 
current  cell  surface.  The  test  segments  are  x-directed  for  1  <  i  <  N2x 
and  y-directed  for  N2x  <  i  <  N2.  The  dot  product  of  the  line  and 
surface  integrals  is  then 
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JC2iffl-J5j]  &,<*«>  •Kt-.j)*' 


Jc2xidXJs3*J^”(i2ij)  T'(l3>)ds/  '  1~i-N 2* 
Jc2yid7  Jsaxj  G™(RU >  I’*(73J)  ds>  '  1  ~ N 2 

Ljx\s3yiGhM  T,(^ds'  : 

JC2yi  dy  JS3yj  Ty(73J)  ds'  ;  N2x<  i  <tf2 


(4-3) 


4.1.1  Impedance  Sub-Matrix  Gaussian  Quadrature  Equations .  Using 
five-node  Gaussian  quadrature,  and  integration  over  the  interval 
x  <  [-1,1]  can  be  approximated  by  the  summation 


d*  =  £  Wn^(xn> 


n-  1 


(4-4) 


where 

wn  -  the  weight  assigned  to  the  nth  Gaussian  quadrature  node  in 
the  interval  [-1,1] 

x„  -  the  location  of  the  nth  Gaussian  quadrature  node  in  the 
interval  [-1,1] 

To  complete  the  triple  integrations  of  equation  (4-1),  the  single 
variable  quadrature  is  extended  to  incorporate  three  variables.  Then, 
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for  example,  the  triple  integrations  of  equation  (4-4)  for  the  index 
ranges  1  <  j  <  N3x  and  1  <  i  i  can  be  approximated  by  the  triple 
summation 


“  £  E  E  WlW«.WnGA"(£ii.£,>llBn)2’x(7jB„) 
1-1  m- 1  n- 1 


(4-5) 


where 

wi.  wm>  wn  “  the  Gaussian  quadrature  weights  as  determined  by  the 
value  of  their  subscript  variables 

1  -  the  quadrature  node  index  for  the  observer  current  cell's  test 
segment 

m,  n  -  the  quadrature  node  indices  for  the  source  current  cell's 
x-coordinate  and  y  coordinate,  respectively 

The  position  vectors  £il  and  £’Jnm  in  equation  (4-5)  now  require 
clarification.  Since  the  Gaussian  quadrature  nodes  and  weights  are 
generally  available  for  the  interval  (-1,1],  a  coordinate  transformation 
must  be  made  with  regard  to  the  integration  variables  of  equation  (4-5) . 
The  appropriate  transformation  is  given  by  [4:2-83] 


r2f<t)dc  -  A2'Ai  r1  f 

A2  +  ■ij  +  Xn(^2“^l) 

J>i  y  ’  2 J-i 

2 

(4-6) 
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On  the  x-directed  observer  current  cells,  the  line  integration 
along  the  cell's  test  segments  means  A2  —  —  a2/ 2  and  A2  -  xA  +  a2/ 2. 

Then,  according  to  equation  (4-6),  the  scale  factor  of  the  coordinate 
transformation  is  a2/2,  and  the  position  vector  is 


y 7i ;  lzi^N 


2z 


(4-7) 


where 

Xi  -  the  1th  Gaussian  quadrature  node  in  the  interval  [-1,1] 

Xj  and  yi  -  the  x-  and  y-coordinates  of  the  position  vector  £XI 

Along  a  y-directed  observer  test  segment,  Ax  -  yt  -  b2/ 2  and 
A2  -  yA  +  b2/2.  The  coordinate  transformation  generates  a  scale  factor 
of  bz/2,  and  the  position  vector  £il  is 


*ii  ~ 


xx< 


Yi 


y  ib 


■  )  ;  ^2x<  i  -^2 


(4-8) 


where 

y1  -  the  1th  Gaussian  quadrature  node  in  the  interval  [-1,1] 
xA  and  y*  -  the  x-  and  y-coordinates  of  the  position  vector  £yi 

The  surface  integrations  require  transformations  in  both  the  x- 
and  y-coordinates.  Over  an  x-directed  current  cell  of  interface  3b,  the 
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intervals  for  the  surface  integrations  are  x  e  [x3  -  a3,  x3  +  a3]  and 
y  €  [y3  -  b3/ 2,  y^  +  b3/2],  Thus,  the  coordinate  transformation  scale 
factor  is  a3b3/2,  and  the  position  vector  £.' is 


Jmn 


-  x  (*j  +  *na3>  +  y  bj- 


(4-9) 


where 

Xn,  -  the  mth  quadrature  node  location  in  the  interval  [-1,1] 
yn  -  the  nth  quadrature  node  location  in  the  intarval  [-1,1] 

Over  the  surface  of  a  y-directed  cell  x  e  [Xj  -  a3/2  ,  x3  +  a3/ 2] 
and  y  t  [y^  -  b3,  y3  +  b3]  .  The  scale  factor  is  again  a3b3/2,  and  3mn 

is  given  by 


*  J. 


-  x 


+  y (Yi  +  ymb3) ; 3 


(4-10) 


With  these  definitions  for  the  position  vectors,  the  variable  7^ 
becomes 


;  lsj<w3x 

y-(£'Jnin ;  n3x<j<n3 


(4-11) 
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Substituting  the  appropriate  coordinate  transformation  scale 
factors  and  summation  approximations  for  the  integrations  of  (4-1) ,  the 
numerical  integration  equations  for  the  impedance  sub-matrix  becomes 


-b23' 
£ij  . 


n2  n3 

£  £1 

i  - 1  k-1 
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where  the  summation  index  j  of  equations  (3-24),  (3-27),  (3-29),  and 
(3-31)  has  been  replaced  with  the  index  /c  to  prevent  confusion  with  the 
use  of  j  to  indicate  the  7-1.  The  other  variables  of  equation  (4-12) 
are  defined  as  follows: 
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4 .1.2  Inner-Produce  Gaussian  Quadrature  Equations.  The  elements 
of  the  forcing  function  in  the  3-plane  MoM  matrix  equation  are  defined 
in  equation  (3-43).  As  with  the  impedance  sub-matrix,  the  first  step  in 
obtaining  the  Gaussian  quadrature  equations  for  the  inner-products  is 
understanding  the  significance  of  the  dot  product  between  the  dyadic 
Green's  function  and  the  magnetic  surface  current  vector  of  the  aper¬ 
ture.  From  equation  (3-39),  the  magnetic  surface  current  vector  is 
y-directed,  therefore  the  dot-product  of  the  general  dyadic  Green's 
function  G^  and  the  magnetic  surface  current  vector  is 


Gb  > 
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y 


(4—14) 
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On  the  x-directed  current  cells,  the  test  functions,  vi(x) ,  are  x- 
directed,  and  over  the  y-directed  current  cells  they  are  y-directed. 
Thus,  the  dot  product  between  the  line  and  surface  integrals  of  equation 
(3-43)  produces  the  following: 

As  with  the  impedance  elements,  Gaussian  quadrature  is  used  to 
numerically  evaluate  the  integrations  of  equation  (4-15)  using  the  same 
summation  index  assignments  as  for  equations  (3-24),  (3-27),  (3-29)  and 
(3-31).  On  the  observer  current  cells,  the  paths  Cxi  and  Cyi  are  such 
that  x  €  [Xj  -  a/2,  Xj^  +  a/2]  and  y  €  fyA  -  5/2,  yA  +  5/2],  respectively. 
Therefore,  the  coordinate  transformation  into  the  interval  [-1,1] 
produces  the  scale  factor  a/2  on  the  x-directed  current  cells  and  5/2  on 
the  y-directed  current  cells.  The  position  vector  is  the  same  as 
defined  in  equations  (4-7)  and  (4-8). 

Since  the  center  of  the  aperture  is  defined  as  the  origin  of  the 
antenna's  coordinate  system,  the  antenna’s  aperture  surface  is  defined 
by  x  <  [-W/2,  W/2]  and  y  e  [-L/2,  L/2 ] .  With  the  required  coordinate 
transformation,  a  scale  factor  of  WL/4  is  obtained,  and  the  source 
position  vector  is  defined 


L  i Is  ib  <£i  • £/ ) My(y/ }  ds' : 

Jc  ;Nx<i<w 


(4-15) 
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(4-16) 


where 

x„  and  yB  -  the  mth  and  nth  Gaussian  quadrature  nodes  in  the 
interval  [-1,1] 

Using  these  definitions,  the  numerical  approximations  for  the 
inner-products  of  equation  (3-42)  are 
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4.2  Software  Overview 

The  MoM  solution  software  is  designed  around  evaluating  the 
summations  of  equations  (4-12)  and  (4-17).  The  operations  of  equation 
(4-17)  are  very  similar  to  those  of  equation  (4-12),  therefore  identify¬ 
ing  the  requirements  for  evaluating  equation  (4-12)  will  also  serve  to 
identify  those  for  equation  (4-17). 
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Equation  (4-12)  must  be  evaluated  to  form  each  of  the  four 
impedance  sub-matrices  that  form  the  MoM  impedance  matrix.  Calculating 
the  elements  of  each  impedance  sub-matrix  requires  the  following  steps: 

(1)  Establish  the  index  of  the  observer  current  cell,  this  is  the 
value  of  summation  index  i  of  equation  (4-12)  and  the  row  index 
for  the  element's  location  within  the  impedance  sub-matrix. 

(2)  For  each  observer  index,  complete  the  triple  summation  over 
indices  2,  m,  and  n  of  (4-12)  for  each  current  cell  on  the  source 
plane.  The  source  current  cell  is  determined  by  summation  index  j 
of  equation  (4-12).  The  index  j  also  determines  the  column  index 
of  the  element's  location  within  the  impedance  sub-matrix. 

Completing  the  triple  summation  over  indices  1,  m,  and  n  in 
equation  (4-12)  requires  the  following  steps: 

(1)  Establish  the  value  of  index  1  which  controls  the  progression 
over  the  quadrature  nodes  of  the  observer  current  cell  test 
segment . 

(2)  For  each  node  on  the  test  segment,  complete  the  summation 
over  indices  m  and  n.  For  each  of  the  possible  combinations  of 
indices  1,  a,  and  n,  |£ilp  -  £'^1  must  be  calculated  according  to 
equations  (4-7)  or  (4-8)  and  (4-9)  or  (4-10)  and  the  surface 
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current  or  surface  charge  density  basis  function  must  be  calculat 
ed.  Based  on  the  value  of  |£U  -  £'Jinn|,  the  appropriate  Green's 
function  must  be  evaluated  using  the  polynomial  model. 

The  difference  between  completing  the  operations  of  equation 
(4-17)  and  (4-12)  is  that  for  equation  (4-17)  the  source  plane  is  now 
the  complete  aperture,  and  the  piecewise  sinusoidal  magnetic  current 
distribution  function  is  used  instead  of  the  rooftop  basis  function. 

Completing  the  antenna  solution  involves  three  basic  functions. 
First,  antenna  parameters  -  operating  frequency,  dielectric  permit¬ 
tivities  and  permeabilities,  element  dimensions  -  and  the  Green's 
functions'  polynomial  coefficients  must  be  read  from  data  files.  Next, 
data  for  calculating  the  integration  node  position  vectors  must  be 
calculated.  The  MoM  matrix  elements  are  then  calculated  and  the  MoM 
matrix  equation  is  solved  to  obtain  the  expansion  coefficients  for  the 
surface  currents  and  surface  charge  densities.  Finally,  the  surface 
currents  are  calculated  using  the  expansion  coefficients. 

The  MoM  software  consists  of  a  main  program  and  13  subroutines. 
The  function  of  the  main  program  is  to  call  subroutines  in  the  appropri 
ate  order  and  set  operation  flags  and  variables  for  the  subroutines. 
Since  the  main  program  does  not  complete  any  major  portion  of  the 
solution,  it  will  not  be  described  in  detail.  Instead,  any  inputs  to 
subroutines  provided  by  the  main  program  will  be  described  as  each 
subroutine  is  described.  The  following  sections  detail  the  subroutines 
that  were  developed  to  complete  the  three  basic  functions  needed  to 
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complete  the  MoM  solution.  Complete  listings  of  the  MoM  software  are 
contained  in  the  Computer  Program  Supplement  for  this  thesis. 


4 .3  Antenna  Data  Subroutines 

Three  subroutines  are  used  to  either  read  data  for  the  solution  or 
calculate  data  needed  in  the  solution.  These  subroutines  are  the  first 
executed  by  the  main  program. 

4 .3.1  Subroutine  INPUT.  The  data  input  is  accomplished  by 
subroutine  INPUT.  Subroutine  INPUT  is  used  to  read  data  from  two 
separate  data  files  and  then  perform  several  short  calculations  with  the 
data.  The  first  data  file  accessed  is  file  GEO. DAT.  This  file  contains 
seven  lines  of  antenna  data.  The  expected  entries  and  read  formats  for 
the  data  file  are  as  follows: 

Line  1)  The  relative  permeabilities  of  the  dielectric 
substrates  between  interfaces  lb  and  2b,  named  MUR1B,  be¬ 
tween  interfaces  2b  and  3b,  named  MUR2B,  and  above  interface 
3b,  named  MUR3B.  All  values  on  the  same  line  in  the  same 
order.  The  read  format  for  these  values  is  five  columns  with 
three  decimal  places. 

Line  2)  The  relative  permittivities  of  the  dielectric 
s  'bstrates  between  interfaces  la  and  2a,  named  EPSRlA,  be¬ 
tween  interfaces  lb  and  2b,  named  EPSR1B,  between  interfaces 
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2b  and  3b,  named  EPSR2B,  and  above  interface  3b,  named 
EPSR3B.  All  values  on  the  same  line  with  the  value  for 
EPSR1A  first.  The  read  format  for  these  values  is  five 
columns  with  three  decimal  places.  The  variables  EPSR1A  and 
EPRSR1B  are  then  used  to  calculate  the  effective  wavenumber, 
variable  KE,  for  the  aperture's  magnetic  field  approxima¬ 
tion. 

Line  3)  The  operating  frequency  of  the  antenna  in  GHz.  The 
read  format  for  the  frequency  is  five  columns  with  two  deci¬ 
mal  places,  thus  limiting  the  maximum  operating  frequency  to 
99.99  GHz.  However,  the  read  format  can  be  changed  to 
accommodate  any  higher  frequencies,  there  are  not  any  soft¬ 
ware  features  that  would  prevent  raising  the  maximum  possi¬ 
ble  frequency.  The  frequency  is  converted  to  radians/second 
and  assigned  the  variable  name  OMEGA. 

Line  4)  The  number  of  integration  nodes  to  be  used  for  all 
line  and  surface  integrations.  The  value  is  read  using  a 
one  column  integer  format  and  is  assigned  the  variable  name 
NUMNODES .  The  software  is  written  for  four  or  six  node 
Gaussian  quadrature. 

Lines  5-7)  These  lines  contain  the  geometry  data  of  the 
antenna.  Line  5  contains  data  for  the  aperture,  Line  6 
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contains  data  for  the  patch  of  interface  2b,  and  Line  7 
contains  data  for  the  patch  of  interface  3b.  Each  line 
contains  the  x-dimension  of  the  aperture  or  patch,  the  y- 
dimension  of  the  aperture  or  patch,  the  height  of  the  inter¬ 
face  above  the  aperture,  and  then  the  x-  and  y-distances 
between  the  patch's  center  relative  to  the  aperture's  cen¬ 
ter.  The  origin  of  the  antenna  is  set  at  the  aperture  cen¬ 
ter,  so  the  height  and  offset  entries  on  Line  4  will  all  be 
zero.  All  dimension  entries  in  the  data  file  are  in  units 
of  centimeters.  After  the  dimension  data  is  read,  it  is 
converted  to  units  of  meters  and  placed  in  the  one-dimen¬ 
sional  arrays  XDIMEN ,  YDIMEN ,  ZPLANE ,  XOFFSET ,  and  YOFFSET . 

The  array  indices  of  these  arrays  designate  to  which  element 
of  the  antenna  the  dimension  applies,  with  one  corresponding 
to  the  aperture,  two  corresponding  to  the  patch  on  interface 
2b,  and  three  corresponding  to  the  patch  on  interface  3b. 

All  dimensional  data  is  read  with  a  format  of  six  columns 
and  three  decimal  places. 

The  second  data  file  accessed  is  file  GREENFS.DAT.  This  file 
contains  the  data  for  the  polynomial  representations  of  the  antenna's 
Green's  functions.  The  coefficients  are  stored  in  3  three-dimensional 
arrays  -  GBA,  IB,  and  GBQ.  The  details  of  the  polynomial  models  for  the 
various  Green's  functions  are  given  in  Appendix  A. 
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Array  GBA  contains  the  polynomial  coefficients  for  the  three 
Green's  functions  that  have  subscripts  Ast,  where  s  and  t  t  [2,3].  The 
first  index  of  an  element  in  array  GBA  designates  which  Green's  function 
the  coefficients  belong  to  -  1  **  G.?->.  2  ♦♦  GA33,  and  3  **  G^  or  G^.  For 
the  A22  Green's  functions,  the  function  with  superscript  bxx  equals  the 
function  with  superscript  byy,  and  the  functions  with  superscript  bxy 
and  byx  are  zero.  Thus,  the  coefficients  in  array  GBA  with  a  first 
index  equal  to  one  represents  all  non-zero  Green's  functions  with 
subscript  A22 .  Similarly,  the  Green's  functions  subscripted  A33  have 
the  bxx  function  equal  to  byy,  with  byx  and  bxy  equal  to  zero.  The 
third  row  of  array  GBA  contains  the  coefficients  of  the  Green's  G^  or 
Ga32 .  For  these  Green's  functions,  the  functions  with  superscript  bxx 
are  equal  to  those  with  superscript  byy,  and  the  functions  with  super¬ 
script  bxy  or  byx  are  zero.  Additionally,  the  Green's  functions  with 
subscript  A23  are  equal  to  those  with  subscript  A32. 

The  Green's  functions  that  have  a  superscript  of  E21  or  E31  are 
defined  by  two  or  more  Sommerfeld  integrals,  and  array  IB  contains  the 
coefficients  for  the  polynomials  representing  these  integrals.  The 
first  index  of  array  IB  determines  which  of  the  six  integral  par^s  the 
coefficients  represent.  Array  elements  with  a  first  index  of  1 ,  3,  and 
5  correspond  to  E21  Green’s  function,  and  the  elements  with  a  first 
index  of  2,  4,  and  6  correspond  to  the  E31  Green's  function. 

Array  GBQ  contains  coefficients  for  the  Green's  functions  needed 
to  calculate  the  electric  fields  created  by  the  surface  charge  density. 


4-17 


The  first  index  of  the  array  corresponds  to  the  subscript  of  the  Green's 
function  —  1  ♦*  q22 ,  2  **  q33 ,  3  «*  q23 ,  and  4  **  q32. 

The  polynomial  models  for  the  GA  and  GE  Green's  functions  all  have 
the  same  form  -  real  component  represented  by  a  polynomial  whose  form  is 
determined  by  the  value  of  R,  and  an  imaginary  component  represented  by 
the  same  polynomial  over  all  values  of  R.  For  these  functions,  a 
second  index  equal  to  1  indicates  coefficients  for  the  polynomial 
modelling  the  real  component  for  R  <  .01,  and  a  second  index  equal  to  2 
indicates  coefficients  for  the  polynomial  modelling  the  real  component 
when  R  >  0.01.  The  elements  with  a  second  coefficient  equal  to  3 
contain  the  coefficients  for  the  polynomial  modelling  the  imaginary 
component  of  the  Green's  functions.  The  third  dimension  of  the  array 
elements  in  arrays  GBA  and  IB  indicates  the  coefficient  subscript  as 
detailed  in  Appendix  A.  The  value  of  the  third  index  is  equal  to  the 
coefficient  subscript  plus  one. 

The  Gq  Green's  functions  are  modeled  with  one  polynomial  which  is 
valid  over  all  values  of  R  representing  the  real  component,  and  differ¬ 
ent  polynomials  valid  over  two  different  intervals  of  R  representing  the 
imaginary  component.  The  second  index  of  array  GBQ  has  the  same 
function  as  the  second  index  of  arrays  GBA  and  IB,  except  an  index  value 
of  one  designates  the  polynomial  for  the  real  component,  a  second  index 
of  two  designates  the  polynomial  data  for  the  interval  R  <  0.01,  and  an 
index  of  3  indicates  data  for  the  interval  R  >  0.01. 

The  coefficients  are  all  real  values  and  are  first  read  into  array 
GBA,  then  IB,  and  GBQ.  Each  Green's  function  has  a  polynomial  model  for 


4-18 


its  real  component  of  imaginary  component  valid  only  for  specific  ranges 
of  R,  and  the  first  line  defines  these  ranges.  The  next  63  lines 
contain  the  data  for  arrays  GBA  and  IB.  The  data  is  arranged  such  that 
the  coefficients  of  the  polynomial  of  the  function's  real  component  for 
R  <  0.01  are  on  the  first  two  lines;  five  then  four  coefficients  per 
line.  Next  come  the  coefficients  for  the  R  >  0.01  real  component 
polynomial;  two  lines  with  five  coefficients  per  line  then  a  single 
coefficient.  Finally,  the  coefficients  for  the  imaginary  component  are 
on  lines  six  and  seven;  five  then  two  coefficients  per  line. 

The  last  28  lines  of  the  GREENFS.DAT  file  contain  the  coefficients 
for  the  Gq  Green's  functions;  seven  lines  per  Green's  function.  The 
first  two  lines  contain  the  coefficients  for  the  real  component's 
polynomial;  five  then  two  coefficients  per  line.  Lines  three  and  four 
contain  the  coefficients  for  the  interval  R  <  0.01  polynomial  of  the 
imaginary  component;  five  then  four  coefficients  per  line.  Lines  five 
through  seven  contain  the  coefficients  of  the  imaginary  component's 
polynomials  for  R  >  0.01;  with  five  coefficients  per  line  on  lines  five 
and  six  and  one  coefficient  on  line  seven.  The  coefficients  are  read 
using  the  format  1PE16 . 7 . 

4 .3.1  Subroutine  SEGMENT.  The  first  data  calculation  subroutine 
is  named  SEGMENT.  This  subroutine  calculates  the  data  needed  to  deter¬ 
mine  the  value  of  R  for  the  surface  and  line  integration  nodes  on  all 
planes . 

As  discussed  in  the  details  of  the  MoM  solution,  each  of  the 
patches  are  divided  into  Nx  x-directed  current  cells  and  Ny  y-directed 
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current  cells.  Each  charge  cell  has  an  x-dimension  of  a  and  a  y-dimen- 
sion  of  b.  Following  Mosig's  work  [17:320],  the  position  vector  data  is 
calculated  using  a  -  0.05A  and  b  such  that  there  are  five  y-directed 
current  cells  for  each  patch.  To  calculate  the  value  of  a,  the  subrou¬ 
tine  assumes  the  x-dimension  of  the  patch  on  interface  3b  is  approxi¬ 
mately  A/2,  thus  a  is  equal  to  one-tenth  of  the  patch's  x-dimension. 
Using  the  a  value  for  interface  3b,  the  patch  on  interface  2b  is  divided 
so  that  its  a  value  is  less  than  or  equal  to  that  of  interface  2b. 

Having  established  the  dimensions  of  the  charge  cells,  the  current  cell 
dimensions  are  known,  and  position  vector  data  can  be  found. 

To  calculate  and  store  the  position  /ector  data,  each  patch's 
current  cells  are  modeled  as  two-dimensional  arrays.  For  example, 

Figure  4-1  illustrates  the  case  of  five  x-directed  charge  cells  and  five 
y-directed  charge  cells  across  the  patch.  In  this  case,  a  4x5  array  of 
x-directed  current  cells  and  a  x5  array  of  y-directed  current  cells  is 
obtained.  The  first  x-directed  current  cell,  in  the  upper  left  corner 
of  the  patch,  corresponds  to  element  (1,1)  of  the  current  cell  array  and 
the  twentieth  x-directed  current  cell,  in  the  lower  right  corner,  corre¬ 
sponds  to  element  (4,5). 

Since  the  current  cell  position  vector  locates  the  cell's  center, 
all  current  cells  in  the  same  column  of  the  array  have  equal  x-compo- 
nents  in  their  position  vectors,  and  all  current  cells  in  i/ue  s*me  row 
of  the  array  have  equal  y-components .  With  the  origin  at  the  lower  left 
of  the  patch,  the  x-  and  y-components  of  the  position  vector  for  any 
current  cell  are 
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First  x-dlrected  currsnt  call 
corresponds  to  srrsy  slement  (1.1) 
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Figure  4-1  Array  Representation  of  Current  Cells  on  the  Antenna  Patch- 
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where 


-  the  array  row  index  of  the  current  cell 

-  the  array  column  index  of  the  current  cell 


-  the  y-dimension  of  the  patch 


On  the  patch  of  Figure  4-1,  the  y-directed  current  cells  form  a 
4x5  array,  with  the  first  y-directed  current  cell  corresponding  to  ele¬ 
ment  (1,1)  being  located  in  the  lower,  left  corner  of  the  patch.  The 
elements  in  the  y-directed  current  cell  array  are  numbered  such  that  a 
larger  column  index  indicates  an  increase  in  the  y-direction,  and  a 
larger  row  index  indicates  an  increase  in  the  x-direction.  With  this 
convention,  the  twentieth  y-directed  current  cell,  corresponding  to 
array  element  (4,5),  is  located  in  the  upper  right  corner  of  the  patch. 
In  this  instance,  the  position  vectors  of  the  current  cells  in  the  same 
array  row  will  have  equal  x-components ,  and  the  position  vectors  of 
current  cells  in  the  same  column  will  have  equal  y-components .  The 
position  vector  components  for  the  y-directed  current  cells  are  given  by 
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To  translate  the  origin  from  the  lower  left  corner  to  the  center 
of  the  patch,  one  half  the  x-dimension  of  the  patch  is  subtracted  from 
each  x-component,  and  one  half  the  y-dimension  of  the  patch  is  subtract¬ 
ed  from  the  y-component.  To  compensate  for  the  possible  offset  between 
the  patch  and  aperture  centers,  the  patch's  center  offset  values  are 
added  to  the  appropriate  coordinates. 

Using  this  array  model,  subroutine  SEGMENT  calculates  the  position 
vector  x-  and  y-components  for  the  rows  and  columns  of  current  cells 
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storing  them  in  4  two-dimensional  arrays.  The  data  arrays  XXCOORD  and 
XYCOORD  contain  the  position  vector  components  for  the  columns  and  rows, 
respectively,  of  x-directed  current  cells  on  the  patch,  and  data  arrays 
YXCOORD  and  YYCOORD  contain  the  same  information  for  the  rows  and 
columns  of  the  y-directed  cells.  The  first  index  of  these  data  arrays 
designates  the  array  row  or  column  in  which  the  cell  is  located,  and  the 
second  index  identifies  on  which  plane  the  cell  is  located. 

The  surface  integration  node  position  vector  components  are 
calculated  using  similar  array  models,  except  the  array  is  contained 
within  an  x-directed  or  y-directed  current  cell.  The  node  array  is 
square,  with  dimensions  being  equ'l  to  the  number  of  Gaussian  quadrature 
nodes ,  variable  NUMNODES .  The  nodes  are  arranged  in  the  array  such  that 
element  (1,1)  corresponds  to  the  node  with  the  most  negative  x-coordi- 
nate  and  most  positive  y-coordinate ,  and  element  (NUMNODES .NUMNODES)  is 
the  node  with  the  most  positive  x-coordinate  and  the  most  negative  y- 
coordinate.  With  the  origin  located  at  the  current  cell’s  center,  on  x- 
directed  current  cells  the  coordinates  for  the  integration  node  array 
rows  and  columns  are 


7  an 


■yab 


(4-20) 


where 
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m  and  n  -  the  summation  indices  as  used  in  equations  (4-12)  and 
(4-17),  m  corresponds  to  the  node  array's  row  index,  n  corresponds 
to  the  node  array's  column  index 

x„  -  the  value  of  the  nth  integration  node  on  x  e  [-1,1] 
ym  -  the  value  of  the  mth  integration  node  on  y  t  [-1,1] 

For  y-directed  current  cells,  the  coordinates  of  the  integration 
node  array  rows  and  columns  are 
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The  position  data  of  the  surface  integration  nodes  is  stored  in  4 
two-dimensional  arrays.  For  nodes  on  an  x-directed  current  cell,  the 
data  array  XCOL  contains  the  x-coordinates  for  the  node  array  columns , 
and  the  data  array  XROW  contains  the  y— coordinates  for  the  rows .  In  a 
y-directed  current  cell,  the  x-coordinates  of  the  rows  in  the  node  array 
are  contained  in  the  data  array  YROW,  and  the  y-coordinates  of  the  node 
array  columns  are  contained  in  data  array  YCOL.  For  all  node  coordinate 
data  arrays,  the  first  index  designates  the  quadrature  node  with  which 
the  coordinate  matches,  and  the  second  index  designates  on  which  plane 
the  node  is  located. 

The  line  integration  nodes  are  modeled  as  a  one-dimensional  array 
along  the  test  segment  of  the  observer  current  cell.  For  an  x-directed 
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current  cell,  the  test  segment  is  centered  in  the  cell  with  respect  to 
the  y-direction,  and  the  segment  covers  the  interval  x  e  [-a/2,  a/2]. 
Since  the  test  segment  is  centered  with  respect  to  the  y-direction,  the 
y-coordinate  for  all  nodes  along  the  test  segment  is  zero,  and  the  x- 
coordinate  is 


*1  " 


nLa 

~r 


(4-22) 


where 

nx  -  the  value  of  the  1th  integration  node  in  the  interval  [-1,1] 

1  -  the  summation  index  as  used  in  (4-12)  and  (4-17) 

On  a  y-directed  current  cell,  the  test  segment  is  centered  with 
respect  to  the  x-direction,  thus  the  x-coordinate  of  all  nodes  along  the 
test  segment  is  zero.  Since  the  y-directed  test  segment  covers  the 
interval  y  e  [-b/2,  b/2] ,  the  y-coordinate  of  the  nodes  along  the  test 
segment  are 


y  i  - 
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The  line  integration  node  coordinates  are  stored  in  2  two-dimen¬ 
sional  data  arrays.  The  data  array  SEGX  contains  the  node  coordinates 
for  the  line  integration  nodes  along  the  test  segment  of  an  x-directed 
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current  cell,  and  SEGY  contains  the  coordinates  for  the  nodes  along  the 
test  segment  of  a  y-directed  cells.  The  first  index  of  the  data  array 
elements  identifies  with  which  node  the  coordinate  belongs,  and  the 
second  index  identifies  the  plane  on  which  plane  the  node  is  located. 

Since  the  antenna  aperture  is  not  divided  into  charge  and  current 
cells,  its  position  vector  calculations  are  slightly  different.  Subrou¬ 
tine  SEGMENT  treats  the  aperture  as  a  single  x-directed  current  cell 
centered  on  the  origin,  thus  the  values  of  XXCOORD(l.l)  and  XYCOORD(l.l) 
are  both  zero.  The  coordinates  of  the  aperture's  integration  nodes  are 


y  oin 
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The  aperture  surface  integration  node  position  vector  components 
are  stored  in  the  data  arrays  XROW  and  XCOL  using  the  same  index  assign¬ 
ment  conventions  as  for  the  nodes  on  the  antenna  patches. 

Since  the  aperture  serves  only  as  a  source  plane  in  the  three 
plane  solution,  there  is  no  need  to  calculate  integration  node  position 
vector  data  for  observer  cell  line  integrations. 

4 .3.2  Subroutine  ASYGFS .  Subroutine  ASYGFS  evaluates  the  value 
of  the  1/R  term  in  the  asymptotic  form  of  the  Green's  functions  G^. 

Ga33,  Gq22  >  anc*  Gq33  which  is  given  by  equation  (3-37).  The  subroutine 
first  calculates  the  scaling  factor  for  the  surface  current,  array 
JFACT,  or  surface  charge,  array  QFACT,  Green's  functions  based  on  the 
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plane  of  interest,  then  the  expression  of  (3-37)  is  evaluated.  The 
total  contribution  of  the  1/R  term  is  calculated  by  multiplying  the 
value  of  equation  (3-37)  by  the  appropriate  scale  function.  The 
contribution  of  the  1/R  term  is  stored  in  two,  2-dimensional  arrays. 
Array  JASYVAL  contains  the  contribution  for  Green's  functions  and 
Ga33.G  The  array  QASYVAL  contains  the  contribution  to  Green's  functions 
Gq22  and  Gq33.  For  both  arrays,  the  first  index  indicates  to  which  plane 
the  contribution  belongs  -  2  «*  interface  2b  and  3  ♦*  interface  3b.  The 
second  index  of  the  arrays  indicates  to  which  orientation  of  cell  the 
contribution  belongs  -  1  **  an  x-directed  current  cell  and  2  **  a  y-direc- 
ted  current  cell. 


4.4  Method  of  Moments  Matrix  Equation  Subroutines 

The  MoM  matrix  equation  is  constructed  and  solved  using  the  first 
level  subroutines  EVALSUB,  EVALRHS  and  SOLVNVRS .  Subroutine  EVALSUB 
evaluates  the  elements  of  each  impedance  sub-matrix  and  places  the 
elements  in  the  proper  locations  of  the  MoM  impedance  matrix.  The 
subroutine  EVALRHS  calculates  the  column  matrix  of  inner  products  that 
form  the  forcing  function  of  the  MoM  matrix  equation.  Subroutine 
SOLVNVRS  is  used  to  solve  the  MoM  matrix  equation  via  inversion  of  the 
impedance  matrix.  All  these  subroutines  rely  on  second  or  third  level 
subroutines  to  complete  their  task.  The  descriptions  of  these  lower 
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level  subroutines  and  their  functions  immediately  follow  the  description 
of  the  higher  level  subroutine  that  uses  them. 

4.4.1  Subroutine  EVALSUB.  As  mentioned  above,  EVALSUB  evaluates 
the  impedance  sub-matrices  of  the  MoH  impedance  matrix.  Subroutine 
EVALSUB  only  calculates  the  impedance  matrix  for  the  observer/source 
plane  combination  designated  by  the  main  program,  thus  four  calls  to 
EVALSUB  are  needed  to  complete  the  MoM  impedance  matrix.  As  each 
element  of  the  sub-matrices  are  evaluated,  EVALSUB  stores  the  element  in 
the  MoM  impedance  matrix,  program  variable  Z.  By  inserting  the  elements 
directly  into  the  MoM  impedance  matrix,  program  memory  requirements  are 
reduced  because  the  need  for  a  large  array  to  store  the  impedance  sub¬ 
matrix  is  eliminated,  and  program  execution  time  is  saved  because  the 
time  required  to  transfer  the  sub-matrices  from  temporary  locations  into 
the  Z  matrix  is  eliminated. 

To  complete  it's  assigned  task,  subroutine  EVALSUB  performs  the 
summations  of  equation  (4-12)  for  each  combination  of  the  observer  plane 
current  cell  index,  variable  name  INDEXOBS,  and  the  source  plane  current 
cell  index,  variable  name  INDEXSRC.  The  variable  INDEXOBS  corresponds 
to  the  summation  index  i  of  equation  (4-12),  and  variable  INDEXSRC 
corresponds  to  the  summation  index  j  of  equation  (4-12).  To  increment 
through  the  observer/source  cell  combinations,  EVALSUB  uses  a  two-level, 
nested  do  loop.  The  outer  loop  controls  variable  INDEXOBS,  and  the 
inner  loop  controls  variable  INDEXSRC.  Since  the  observer  plane  current 
cell  index  controls  the  row  index  of  the  element  in  the  impedance  sub¬ 
matrix  and  the  source  cell  index  controls  the  column  index,  these  two  do 


4-28 


loops  actually  determine  which  element  of  the  sub-matrix  is  being 
evaluated. 

Before  the  summations  are  started,  EVALSUB  uses  row  and  column 
indices  of  the  present  impedance  matrix  element,  the  appropriate  scale 
factor,  variable  SCALEC,  for  the  summations 


SCALEC  - 


•  / 

jwaa 

— zr~ 

jubb' 

— 1 r~ 


;  1NDEXSRC INDEXOBS  SNX 
;  INDEXSRC>N'x,  INDEXOBS>Nx 


(4-25) 


where  the  prime  indicates  source  plane  variables. 

Since  the  Green's  functions  with  superscripts  byx  and  bxy  in 
equation  (4-12)  are  zero,  the  triple  summations  of  the  equation  are  only 
calculated  when  the  observer  and  source  cells  are  both  x-directed  or  y- 
directed  cell  and  the  value  of  SCALEC  is  only  needed  for  the  ranges  of 
indices  i  and  J  as  indicated  in  equation  (4-25).  The  triple  summations 
of  equation  (4-12)  are  calculated  in  part  by  EVALSUB  and  in  part  by  a 
second  level  subroutine  named  SURFINT.  Subroutine  EVALSUB  completes  the 
summation  over  the  index  1  of  equation  (4-12),  and  SURFINT  evaluates  all 
summations  for  the  indices  m  and  n  of  equation  (4-12).  The  details  of 
subroutine  SURFINT  are  presented  section  4.4.2. 

To  complete  the  triple  summations  of  equation  (4-12),  EVALSUB  uses 
a  do  loop  that  is  executed  once  for  each  node  on  the  observer  cell  test 
segment.  At  each  node,  SURFINT  is  called  to  calculate  the  double  summa¬ 
tion  over  indices  m  and  n  returning  the  result  in  variable  VALUSURF. 
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Subroutine  EVALSUB  then  sums  the  product  of  VALUSURF,  the  appropriate 
Gaussian  node  weight,  and  SCALEC  for  all  nodes  on  the  test  segment.  The 
sum  for  all  test  segment  nodes  is  assigned  the  variable  name  CELEMENT. 

After  completing  the  evaluation  of  the  surface  current  field 
contributions,  the  fields  produced  by  the  surface  charge  density  are 
calculated  with  two  calls  to  SURFINT.  The  first  call  evaluates  the  sum 
over  indices  m  and  n  with  the  observer  point  located  at  £*xi,  and  the 
second  evaluates  contributions  for  the  observer  point  at  .  The 
variable  LQ  is  used  to  inform  subroutine  SURFINT  that  the  call  is  for 
evaluation  of  a  charge  contribution.  The  flag  LQ  can  have  three  values, 
LQ  -  0  implies  the  call  to  SURFINT  is  for  evaluation  of  surface  current 
field  contributions,  LQ  -  1  implies  the  call  is  for  surface  charge  field 
contributions  at  £+Xj ,  and  LQ  -  2  implies  the  call  is  for  surface  charge 
field  contributions  at  £_~xi .  Subroutine  SURFINT  returns  the  value  of 
the  double  summation  in  variable  VALUSURF,  then  EVALSUB  multiplies  the 
value  by  SCALEQ  -  (j2u>)-1.  The  contributions  of  the  surface  charge 
density  are  placed  in  variable  QELEMENT.  The  final  value  of  the  sub¬ 
matrix  element  is  found  by  adding  variables  CELEMENT  and  QELEMENT. 

The  sub-matrix  element  is  then  placed  in  the  proper  element  loca¬ 
tion  of  Z,  the  MoM  impedance  matrix.  Subroutine  EVALSUB  calculates  the 
Z  matrix  location  for  each  sub-matrix  element  using  the  variables 
ISTRTROW,  ISTRTCOL,  INDEXOBS,  and  INDEXSRC.  The  variables  ISTRTROW  and 
ISTRTCOL  identify  the  MoM  impedance  matrix  row  and  column  in  which  the 
first  element  of  the  current  impedance  sub-matrix  is  located.  Thus,  the 
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MoM  impedance  matrix  row  index  is  (ISTRTROU  -  INDEXOBS  -  1).  The  column 
index  for  the  sub-matrix  element  is  (ISTRTCOL  -  INDEXSRC  -  1). 

4.4 .2  Subroutine  SURFINT.  The  second  level  subroutine  SURFINT  is 
used  by  first  level  subroutines  EVALSUB  and  EVALRHS  to  calculate  the 
variable  VALUSURF.  When  called  by  subroutine  EVALSUB,  SURFINT  performs 
the  double  summations  over  the  indices  m  and  n  in  equation  (4-12).  When 
called  by  EVALRHS,  SURFINT  calculates  the  double  summation  over  indices 
m  and  n  of  equation  (4-17) . 

To  complete  the  double  summations,  SURFINT  uses  a  two-layer, 
nested  do  loop.  The  outer  loop  controls  summation  index  m  in  equations 
(4-12)  and  (4-17)  through  variable  MSURFY.  As  implied  by  its  name, 
MSURFY  controls  the  y-coordinate ,  or  the  node  array  column  index,  of  the 
source  cell  surface  integration  nodes.  The  inner  do  loop  controls  the 
summation  index  n  of  equation  (4-12)  and  (4-17)  through  variable  NSURFX. 
Index  NSURFX  controls  the  x-coordinate ,  or  the  node  arr  row  index,  of 
the  source  cell  surface  integration  nodes.  The  number  of  times  these 
two  do  loops  are  executed  is  determined  by  the  value  of  NUMNODES . 

When  called  by  EVALSUB  to  find  fields  from  surface  currents, 
SURFINT  completes  the  double  summations  over  the  source  cell  surface 
integration  nodes  for  the  present  observer  test  segment  node.  SURFINT 
calls  the  third  level  subroutine  RHOFIND  to  calculate  R  for  the  present 
source  surface  integration  and  observer  test  segment  nodes.  Subroutine 
RHOFIND  uses  the  variables  INDEXOBS,  INDEXSRC,  MSURFY,  and  NSURFX  and 
appropriate  elements  of  the  node  coordinate  data  arrays  to  determine  the 
coordinates  for  the  observer  and  source  plane  integration  nodes.  Details 
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of  RHOFIND  are  presented  in  section  4.4.3.  The  value  of  R  is  returned 
in  variable  RHODIFF.  Subroutine  RHOFIND  also  calculates  the  rooftop 
function,  variable  CAPT  Next,  SURFINT  calls  subroutine  EVALGRN  to  find 
value  of  the  appropriate  Green's  for  the  present  value  of  RHODIFF  which 
is  returned  in  variable  GF.  Details  of  subroutine  EVALGRN  are  given  in 
section  4.4.4. 

Next,  the  product  of  the  Gaussian  node  weights  for  MSURFY  and 
NSURFX  is  computed  and  given  the  variable  name  GAIN.  If  the  call  from 
EVALSUB  was  for  evaluation  of  a  surface  current  contribution,  the 
variable  LQ  will  be  set  to  0  by  EVALSUB,  and  the  product  of  GAIN,  CAPT, 
and  GF  is  summed  for  all  surface  integration  nodes  to  obtain  VALUSURF. 

If  the  call  from  EVALSUB  is  for  the  evaluation  of  a  surface  charge 
contribution,  the  value  of  variable  LQ  will  be  1  or  2  and  VALUSURF  is 
calculated  differently.  In  this  case,  subroutine  RHOFIND  will  calculate 
RHODIFF  and  the  surface  charge  density  basis  function,  variable  PI,  for 
the  observer  being  located  at  the  appropriate  end  of  the  observer 
current  cell  test  segment.  The  Green's  function  for  the  charge  density 
contributions  will  then  be  evaluated,  and  the  product  of  GAIN,  PI,  and 
GF  will  be  summed  to  obtain  VALUSURF. 

When  called  by  EVALRHS ,  SURFINT  completes  the  summations  of 
equation  (4-17)  for  the  aperture  surface  integration  and  the  observer 
test  segment  nodes.  In  this  case,  the  do  loop  indices  MSUR1-".  and  NSURFX 
control  the  row  and  column  indices  of  the  array  of  integration  nodes  on 
the  surface  of  the  aperture.  The  call  to  subroutine  RHOFIND  will  return 
RHODIFF;  the  value  of  the  piecewise  sinusoidal  function  of  equation 
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(3-37),  variable  PWS;  and  the  angle  f,  variable  ZETA,  as  defined  in 
equation  (3-5)  for  y-directed  sources.  The  appropriate  Green's  function 
polynomial  and  GAIN  are  then  evaluated,  and  VALUSURF  is  obtained  by 
summing  the  product  of  GAIN,  PWS,  and  GF  for  all  nodes  on  the  aperture. 

If  the  source  and  observer  are  coplanar  and  in  the  same  current 
cell,  the  contribution  of  the  1/R  term  as  calculated  by  subroutine 
ASYGFS  is  added  to  VALUSURF. 

4.4 .3  Subroutine  RHOFIND.  Subroutine  RHOFIND  is  a  third-level 
subroutine  that  is  called  by  the  second  level  subroutine  SURFINT.  Sub¬ 
routine  RHOFIND  calculates  R  for  the  source  and  observer  integration 
nodes,  the  surface  current  or  charge  density  basis  function  or  the 
piecewise  sinusoidal  function  for  the  source  integration  node,  and  the 
angle  f . 

Subroutine  RHOFIND  first  finds  the  position  coordinates  for  the 
source  integration  node.  The  first  task  is  to  find  the  position  of  the 
current  cell  relative  to  the  plane's  center.  Next,  the  position  of  the 
present  integration  node  relative  to  the  cell's  center  is  found.  The 
integration  node’s  position  relative  to  the  antenna's  origin  is  then 
found  using  the  relative  positions  of  the  cell  and  node. 

Using  the  array  model  to  represent  the  current  cell's  on  the 
planes  of  the  antenna  means  finding  the  position  coordinates  of  any  cell 
is  accomplished  by  first  determining  the  array  row  and  column  indices 
for  the  cell  and  then  extracting  position  data  from  the  appropriate  data 
arrays,  XXCOORD  and  XYCOORD  or  YXCOORD  and  YYCOORD . 
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The  array  indices  are  found  by  using  variable  INSEXSRC  to  deter¬ 
mine  the  cell's  ordinal  number  within  the  array  of  x-directed  or 
y-directed  current  cells.  Recall,  1  <  INDEXSRC  <  N,  where  N  is  the 
total  number  of  x-  and  y-directed  cells  on  the  source  plane.  The 
current  cells  on  a  plane  are  numbered  such  that  INDEXSRC  -  1  for  the 
first  x-directed  cell  and  INDEXSRC  -  N  for  the  last  y-directed  cell. 
Therefore,  for  an  x-directed  current  cell,  INDEXSRC  is  the  cell's 
ordinal  number. 

Once  the  ordinal  number  is  known,  the  row  and  column  indices  are 
found  in  two  steps.  First,  the  ordinal  number  is  divided  by  the  number 
of  current  cells  per  row  of  the  x-directed  cell  array.  The  whole  number 
portion  of  the  quotient  is  one  greater  than  the  cell's  row  index,  vari¬ 
able  IROW.  For  instance  with  4  cells  in  a  row  and  an  ordinal  number  of 
9,  the  quotient  is  2.25,  thus  IROW  -  3. 

Next  the  column  index,  variable  ICOL,  is  found  by  subtracting  the 
product  of  the  row  index  less  one  and  the  number  of  cells  in  a  row  of 
the  x-directed  cell  array  from  the  cell's  ordinal  number.  Continuing 
the  example,  with  the  row  index  of  3,  ICOL  -  9  -  [4  x  (3-1)]  -  1.  Thus, 
for  the  patch  of  Figure  8,  the  ninth  x-directed  cell  corresponds  to 
element  (3,1)  of  the  current  cell  array.  The  x-coordinate  of  this 
cell's  center  is  contained  in  data  array  element  XXC00RD(3 , JSRCPLAN) , 
and  the  cell  center's  y-coordinate  is  contained  in  XYC00RD(1 .JSRCPLAN) . 

To  obtain  the  ordinal  number  for  a  y-directed  cell,  the  number  of 
x-directed  cells  on  the  plane,  variable  NX ( JSRCPLAN) ,  must  be  subtracted 
from  the  value  of  INDEXSRC.  The  Indices  IROW  and  ICOL  are  then  calcu- 
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lated  in  the  same  steps  as  those  for  an  x-directed  cell.  The  x-coordi- 
nate  of  the  cell's  center  is  YXCOORD ( IROW , JSRCPLAN)  and  the  y-coordinate 
of  the  cell's  center  is  YYCOORD(ICOL, JSRCPLAN) . 

Since  SURFINT  progresses  through  the  array  of  integration  nodes  on 
the  current  cell  using  a  row  and  column  index,  no  manipulations  of  the 
integration  node  indices  are  needed.  For  all  cells,  the  x-coordinate  of 
the  node's  relative  position  is  XCOL(N.jURFX,  JSRCPLAN)  and  the  y-coord¬ 
inate  is  YROW(MSURFY, JSRCPLAN 

The  cell  and  node  relative  position  coordinates  are  summed  to  find 
the  node's  absolute  position  coordinates  and  the  values  are  stored  in 
the  one  dimensional  array  RHOPRIME.  Element  RHOPRIME(l)  contains  the 
x-coordinate  of  the  node's  location,  and  the  element  RH0PRIME(2) 
contains  the  y-coordinate. 

Having  determined  the  source  node  position,  the  basis  functions  of 
equations  (3-13)  and  (3-18)  and  the  magnetic  current  of  equation  (3-40) 
can  be  calculated.  From  equations  (3-13)  and  (3-18),  the  basis  func¬ 
tions  are  dependent  on  the  value  of  (£'j  -  £,^mn),  which  is  actually  the 
position  of  the  integration  node  relative  to  the  center  of  the  cell. 
Thus,  the  surface  current  basis  function,  CAPT,  and  the  surface  charge 
density  basis  function,  PI,  are  calculated  as  follows 
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CAPT 


r.  XCOL (NSURFX , JSRCPLAN) 
A (JSRCPLAN) 

.  XCOL (NSURFX, JSRCPLAN) 
A (JSRCPLAN) 


;  XCOL (NSURFX , JSRCPLAN )  <  0 


INDEXSRC 


; XCOL (NSURFX, JSRCPLAN)  >  0 


<  NX (JSRCPLAN) 


1 ; Otherwise 
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,  .  YROW (MSURFY , JSRCPLAN) 
B (JSRCPLAN) 

,  YROW ( MSURFY , JSRCPLAN ) 
B (JSRCPLAN) 


; YROW (MSURFY, JSRCPLAN) <  0 


INDEXSRC 


;  YROW  (MSURFY,  JSRCPLAN)  >0; 


<  NX (JSRCPLAN) 


1 ;  Otherwise 


PI 


f-l  l  XCOLtNSURFX, JSRCPLAN)  <0 
1  ;  XCOL (NSURFX, JSRCPLAN) >0 


INDEXSRC <  NX (JSRCPLAN) 


0 ;  Otherwise 


-1 

1 

.  0 


YROW  (MSURFY,  JSRCPLAN)  <  O' 
YROW (MSURFY , JSRCPLAN) >  0 


INDEXSRC <  NX (JSRCPLAN) 


Otherwise 
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If  the  source  plane  is  the  antenna  aperture,  then  the  piecewise 
sinusoidal  magnetic  current  function,  variable  PWS,  is  evaluated  by 
RHOFIND.  Since  the  aperture  is  treated  as  a  single  x-directed  current 
cell,  the  integration  node's  y-coordinate  on  the  aperture,  variable 
XROW ( MSURFY, 1) ,  determines  the  value  of  PWS.  Subroutine  SEGMENT 
segments  the  aperture  such  that  A(l)  «  W/2  and  B(l)  -  L,  therefore 


PWS 


sin  ^C£« 


3  (p  ~  I  MSURFY ( I ROW , JSRCPLAN)  | 


2.A(l).«in*g^<l> 


(4-28) 


where  *  denotes  scalar  multiplication. 
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The  observer  test  segment  node  location  is  found  in  the  same 
way  as  the  source  integration  node  location,  except  the  ordinal  number 
of  the  cell  is  found  based  on  the  value  of  INDEXOBS,  and  the  data 
arrays'  second  indices  are  now  IOBSPLAN. 

When  the  field  contributions  from  surface  currents  are  being 
evaluated,  the  observer  cell  integration  node  index  is  controlled  by 
variable  LINENODE.  Since  the  test  segments  are  centered  in  the  current 
cell,  only  one  coordinate  is  needed  to  locate  the  node  relative  to  the 
current  cell's  center.  On  an  x-directed  cell,  the  line  integration  is 
along  an  x-directed  contour,  and  the  coordinates  of  the  node  are 


RH0(1)  -  XXCOORD(ICOL, IOBSPLAN)  +  SEGX(LINENODE , IOBSPLAN) 
RH0(2)  -  XYCOORD(IROW, IOBSPLAN) 


On  a  y-directed  cell,  the  line  integration  is  along  a  y-directed 
contour,  and  the  coordinates  of  the  node  are 


RH0(1)  -  XXCOORD(IROW, IOBSPLAN) 

RH0(2)  -  YYCOORD(ICOL, IOBSPLAN)  + SEGY(LINENODE , IOBSPLAN) 


If  the  contributions  from  the  surface  charge  density  is  being 
evaluated,  the  observer  location  is  at  either  end  of  the  cell's  test 
segment.  The  value  of  variable  LQ  informs  RHOFIND  at  which  end  of  the 
test  segment  the  observer  is  located  as  discussed  in  section  4. 4.1.  The 
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location  of  the  test  segment's  ends  relative  to  the  current  cell's 
center  is  given  in  equation  (11).  Therefore,  on  an  x-directed  cell  the 
x-coordinate  of  the  test  segment  ends  are  found  as  follows: 


RHO(l)  - 


XXCOORD (ICOL, IOBSPLAN)  ♦  3  (IOBSPLAN) 

2 

XXCOORD (ICOL, IOBSPLAN)  -  ft (IOBSPLAN) 


LQ-  1 


LQ-2 


(4-31) 


and  the  y-coordinate  for  the  end  of  the  test  segment  is  equal  to  that  of 
the  current  cell's  center. 

For  a  y-directed  current  cell,  the  x-coordinate  of  the  test 
segment's  ends  is  equal  to  that  of  the  current  cell,  and  the  y-coordi¬ 
nate  for  the  ends  is  given  by 


YYCOORD (ICOL, IOBSPLAN) 


RHO  (2)  - 


YYCOORD (ICOL, IOBSPLAN) 


B (IOBSPLAN)  .  LQ  _  x 

B (IOBSPLAN)  ,  „ 

- -  t  *-y  • 


(4-32) 


Having  calculated  the  x-  and  y-coordinates  for  the  positions  of 
the  observer  and  source  integration  nodes,  the  value  of  variable 
RHODIFF,  can  be  found 


i 

RHODIFF  -  [(RHO  (1)  -  RHOPRIME <1 )  )2  *  (RHO  (2)  -  RHOPRIME  (2)  )J]  J 


(4-33) 
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If  the  aperture  is  the  source  plane,  JSRCPLAN  -  1  and  the  variable 
ZETA  is  calculated  in  accordance  for  equation  (3-5)  for  y-directed 
sources.  Since  the  x-  and  y-coordinates  of  the  position  vectors  are 
known,  the  value  of  ^  and  4>'  are  the  inverse  tangent  of  the  y-coordinate 
divided  by  the  x-coordinate . 

4.4.4  Subroutine  EVALGRH .  Subroutine  EVALGRN  is  a  third  level 
subroutine  called  by  SURFINT  to  calculate  the  value  of  the  appropriate 
Green's  function  polynomial  given  the  present  value  of  RHODIFF. 

First,  value  of  RHODIFF  is  compared  to  the  range  information 
obtained  from  data  file  GREENFS  DAT  to  determine  which  coefficients  and 
polynomial  form  are  to  be  used.  Next,  the  subroutine  determines  which 
Green's  function  is  needed  based  on  the  values  of  IOBSPLAN,  JSRCPLAN, 
INDEXOBS ,  NXv IOBSPLAN) ,  and  LQ .  Calculating  the  Green's  function  is 
then  only  a  matter  of  evaluating  the  appropriate  polynomials  to  obtain 
the  components  of  the  Green's  functions’  value.  The  result  of  the 
calculations  is  returned  to  SURFINT  in  variable  GF. 

4.4.4  Subroutine  EVALRHS.  Subroutine  EVALRHS  is  the  first-level 
subroutine  that  calculates  the  elements  in  the  column  matrix  of  inner- 
products  as  defined  in  equation  (4-17).  Subroutine  EVALRHS  is  called 
twice  by  the  main  program.  The  first  call  evaluates  the  inner  products 
for  the  patch  of  interface  2b  and  the  second  evaluates  the  inner- 
products  for  the  patch  of  interface  3b. 

In  equation  (4-17),  the  matrix  element’s  row  is  determined  by  the 
observer  current  cell  index,  with  one  row  for  each  cell  on  the  observer 
plane.  The  source  plane  for  all  calculations  is  the  antenna  aperture. 
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Since  the  aperture  is  modeled  as  a  single  x-directed  current  cell,  the 
value  of  INDEXSRC  is  always  one. 

Subroutine  EVALRHS  controls  the  variable  INDEXOBS  to  calculate  the 
inner-product  for  each  observer  cell.  At  each  observer  cell,  EVALRHS 
completes  the  summations  of  equation  (4-17)  by  controlling  the  summation 
over  index  2,  variable  LINENODE,  and  calling  subroutine  SURFINT  to  com¬ 
plete  the  summation  over  indices  a  and  n  of  equation  (4-17).  For  calls 
from  EVALRHS,  SURFINT  returns  the  variable  VALUSURF  where 


VALUSURF 


HUMNODES  NUMNODES 


£ 


E 


an)My(y')  ;  1  <INDEXOBS<N2x(4-34) 

)Mw(y')  ;  N,_<  INDEXOBS<N 

on  j  * 


Subroutine  EVALRHS  then  uses  VALUSURF  to  complete  the  calculation 
of  the  inner-product  as  follows: 


RHS (ISTRTR0W+INDEX0BS-1) 


- 


aWL 

T" 

bWL 


NUMNODES 

£  VALUSURF 


l-i 

NUMNODES 


s 


VALUSURF 


(4-35) 


The  result  of  equation  (4-35)  is  then  placed  in  the  appropriate 
location  of  the  inner  product  data  array,  variable  name  RHS.  To  inform 
EVALRHS  of  the  correct  element  index  of  RHS,  the  main  program  sets 
variable  ISTRTROW  before  EVALRHS  is  called  each  time.  Variable  ISTRTROW 
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identifies  in  which  RHS  element  the  inner-product  for  the  first  observer 
cell  must  be  stored.  The  index  for  the  element  of  matrix  RHS  being 
calculated  is  given  by  (ISTRTROW  +  INDEXOBS  -  1). 


4.5  Matrix  Equation  Solution  Subroutines 

Subroutine  SOLVNVRS  solves  for  the  expansion  coefficients  using 
matrix  inversion  and  multiplication.  The  inversion  of  the  MoM  impedance 
matrix  is  completed  by  a  second-level  routine  named  INVERT. 

4 .5.1  Subroutine  SOLVNVRS.  Subroutine  SOLVNVRS  finds  the 
solution  to  the  general  matrix  equation  Ax  -  b,  x  being  the  unknown 
vector,  by  calculating  x  -  A-1b.  The  subroutine  call  for  SOLVNVRS  is 
"CALL  SOLVNVRS (Z,  IDIM,  RHS,  ALPHAS)”  where  Z  is  the  IDIM  x  IDIM  matrix 
equivalent  to  the  matrix  A  in  the  equation  Ax  -  b,  RHS  is  the  known 
forcing  function  of  the  equation  equivalent  to  the  vector  b,  and  ALPHAS 
is  the  unknown  equivalent  to  the  vector  X- 

Subroutine  SOLVNVRS  first  calls  subroutine  INVERT  to  find  the 
inverse  of  Z  which  is  returned  in  the  matrix  ZINV.  Next,  the  solution 
matrix  ALPHAS  is  calculated  by  forming  the  matrix  product  of  ZINV  and 
RHS. 

4 .5.2  Subroutine  INVERT.  Subroutine  INVERT  is  a  third-level 
subroutine  called  by  subroutine  SOLVNVRS  to  calculate  the  inverse  of  the 
matrix  specified  in  the  subroutine  call.  The  subroutine  call  is  'CALL 
INVERT ( BINPUT ,  IDIM,  INVERSE)  where  BINPUT  is  the  matrix  to  be  inverted, 


4-41 


IDIM  is  the  dimension  of  BINPUT,  and  INVERSE  is  the  variable  name 
assigned  to  the  inverse  of  BINPUT.  Subroutine  INVERT  uses  Gaussian 
elimination  to  calculate  the  inverse  matrix. 

Since  analysis  of  the  antenna  requires  all  original  matrices  be 
preserved,  the  first  task  of  INVERSE  is  to  transfer  the  matrix  BINPUT 
into  a  local  matrix  INPUT.  Duplicating  the  BINPUT  matrix  is  necessary 
because  the  reduction  to  triangular  form  will  overwrite  elements  of  the 
original  matrix  with  row  permutation  parameters.  Since  the  inverse 
matrix  will  be  calculated  column  by  column,  the  alterations  to  the 
original  matrix  must  be  tracked  via  permutation  parameters  so  the  same 
alterations  may  be  made  to  each  identity  matrix  column.  More  detail  on 
the  permutation  parameters  is  in  the  subsequent  discussion. 

Once  the  BINPUT  matrix  is  copied,  the  next  step  is  to  transform 
the  INPUT  matrix  into  triangular  form  using  the  standard  row  manipula¬ 
tions  of  Gaussian  elimination.  Partial  pivoting  is  used  to  reduce 
rounding  errors  during  Gaussian  elimination  [12:28].  As  the  rows  of  the 
matrix  are  pivoted,  the  pivot  changes  are  stored  in  the  one-dimension¬ 
al  array  IPIVOT.  For  instance,  if  row  I  was  pivoted  with  row  J  during 
the  elimination  process,  IPIVOT(I)  would  equal  J.  If  the  row  was  not 
pivoted,  the  element  IPIVOT(I)  will  be  I. 

As  each  row  is  changed,  row  permutation  parameters  are  stored  in 
the  INPUT  matrix  elements  that  would  be  zero  with  the  matrix  in  lower 
triangular  form.  For  instance,  if  the  element  INPUT(I,J)  is  reduced  to 
zero  by  multiplying  all  elements  of  row  J  by  the  factor  INPUT(I , J)/INPU- 
T(J,J)  and  then  subtracting  the  result  from  all  elements  of  row  I,  the 
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permutation  parameter  for  INPUT(I,J)  is  -INPUT(I , J)/INPUT(J , J) .  The 
changes  to  every  row  are  recorded  in  this  manner. 

After  the  matrix  is  in  triangular  form,  the  row  permutation 
parameter  contained  in  element  INPUT (K.L)  will  show  what  factor  row  L 
was  multiplied  by  so  that  subtracting  row  L  from  row  K  would  produce  a 
zero  value  for  element  INPUT (K,L). 

With  the  matrix  in  triangular  form,  the  inverse  matrix  is  calcu¬ 
lated  a  column  at  a  time  by  backsolving  the  system  Ax  -  b  where  A  is  the 
triangular  matrix  INPUT,  £  is  the  ith  column  of  the  inverse  matrix,  and 
b  is  the  ith  column  of  the  IDIM  x  IDIM  identity  matrix,  local  variable 
RHS . 

Before  backsolving  for  the  inverse  column  elements,  matrix  RHS 
must  be  altered  in  the  same  way  as  the  original  INPUT  matrix  was  altered 
during  reduction  to  triangular  form.  The  elements  of  RHS  are  trans¬ 
formed  one  at  a  time  by  repeating  these  three  steps  for  each  element. 

1)  First,  as  the  element  RHS(K)  is  altered,  the  entry  for 
IPIVOT(K)  is  checked,  if  IPIVOT(K)  *  K,  then  elements  RHS (K) 
and  RHS ( IPIVOT (K) )  are  exchanged. 

2)  Sum  the  products  of  the  current  entry  in  RHS(K)  and  the 
permutation  parameters  located  on  row  K  of  matrix  INPUT 
which  have  column  indices  less  than  K.  This  summation, 
variable  name  DEL(K) ,  is  the  total  change  for  RHS(K)  result¬ 
ing  when  matrix  INPUT  was  reduced  to  triangular  form. 
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3)  Form  the  transformed  RHS(K)  by  adding  DEL(K)  to  RHS(K). 


After  these  manipulations,  the  column  matrix  RHS  is  in  the  same 
form  as  if  it  had  been  placed  in  an  augmented  matrix  with  INPUT  and 
undergone  the  row  manipulations  while  INPUT  was  being  reduced  to 
triangular  form. 

When  RHS  has  undergone  the  appropriate  changes,  the  system  can  be 
backsolved  to  obtain  the  elements  for  the  column  of  the  inverse  matrix. 
The  inverse  matrix  elements  are  found  beginning  with  element  IDIM  and 
continuing  to  element  1.  The  algorithm  used  to  complete  the  backsolving 
is  the  same  as  that  found  in  Figure  2-1  of  [10],  except  the  variable 
names  have  been  changed  to  match  those  of  subroutine  INVERSE,  and  as 
each  element  of  the  inverse  matrix  is  found,  the  element  is  placed  into 
its  proper  position  in  the  matrix  INVERSE. 


4.6  Data  Output  Subroutines 

All  data  output  from  the  MoM  solution  software  is  handled  by  the 
first-level  subroutine  OUTPUT.  Subroutine  OUTPUT  uses  the  second-level 
subroutine  JDISTR  to  calculate  data  for  the  surface  currant  distribu¬ 
tions  . 

4 .6.1  Subroutine  OUTPUT.  Subroutine  OUTPUT  first  writes  the 
surface  current  expansion  coefficients  to  file  COEFFS.OUT.  Before 
writing  to  file,  OUTPUT  performs  several  calculations  to  determine  the 
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format  to  be  used  for  the  coefficients.  The  coefficients  are  written  in 


two  columns,  with  the  coefficients  for  the  x-directed  then  y-directed 
cells  on  the  patch  of  interface  2b  being  the  first  in  the  file.  After 
the  coefficients  are  written  to  file,  the  subroutine  to  calculate  the 
surface  current  distributions  is  called. 

U.6.2  Subroutine  JDISTR.  Subroutine  JDISTR  calculates  the  value 
of  the  surface  currents  along  the  centers  of  each  row  of  x-directed  and 
y-directed  current  cells  on  both  antenna  patches.  For  each  current 
cell,  the  surface  current  function  is  calculated  at  five  sample  points. 

The  surface  current  data  is  first  calculated  for  the  cells  on  the 
patch  of  interface  2b  then  interface  3b,  with  the  x-directed  cell  data 
being  the  first  calculated  for  both  patches.  There  are  five  samples 
taken  at  equal  intervals  across  each  current  cell.  For  instance,  for  x- 
directed  cells,  the  sample  points  are  spaced  a  distance  of  a/2  with  the 
first  point  being  located  on  the  leftmost  boundary  of  the  left  charge 
cell.  With  these  sample  point  locations  and  the  overlap  of  current 
cells,  the  three  sample  points  in  the  right  charge  cell  of  the  current 
cell  overlap  with  the  first  three  sample  points  in  the  left  charge  cell 
of  the  current  cell  immediately  to  the  right.  Subroutine  JDISTR 
compensates  for  the  overlap  by  summing  the  surface  current  contributions 
for  both  cells  sharing  a  charge  cell. 

The  progression  over  a  row  of  charge  cells  is  controlled  by  a  two 
layer,  nested  do  loop.  The  outer  most  loop  controls  the  cell  within  the 
row  and  the  inner  do  loop  controls  the  index  of  the  sample  point.  At 
each  cell  in  the  row,  the  index  of  the  expansion  coefficient  for  the 
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cell  is  first  calculated,  then  the  surface  current  value  at  the  sample 
point  is  found. 

After  the  surface  current  values  for  the  complete  row  have  been 
calculated,  they  are  normalized  to  units  of  milliamps/square  meter. 

Next,  the  y-coordinate  (x-coordinate)  of  the  x-directed  (y-directed) 
row's  center  is  written  to  file  along  with  the  x-coordinate  (y-coordi¬ 
nate)  and  surface  current  value  at  each  sample  point.  Since  the  surface 
current  is  a  complex  quantity,  the  real  components  are  written  in  file 
J REAL. OUT  and  the  imaginary  components  are  written  to  file  J1MAG.0UT. 
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V.  Results 


This  chapter  presents  the  data  obtained  from  the  MoM  solution 
software  described  in  Chapter  IV.  The  MoM  software  was  used  to  find  the 
solutions  for  two  different  antennas  using  both  four  and  six  node 
Gaussian  quadrature.  The  results  for  each  run  of  the  software  are 
contrasted  and  compared. 


5.1  Physical  Description  of  the  Sample  Antennas 

The  MoM  solution  data  presented  in  this  chapter  was  obtained  using 
the  following  examples  of  aperture-fed,  stacked  patch  antennas.  All 
solutions  are  for  an  operating  frequency  of  3.7  GHz. 

5.1.1  Antenna  Pi.  The  antenna  example  referred  to  as  antenna  *1 
has  resonance  frequencies  of  3.62  Ghz  and  4.05  GHz  [30:48]  and  is 
described  as  follows: 

-  The  thickness  of  each  substrate  is  0.158  cm,  and  all  substrates 

have  tr  «  2.2  and  pz  -  1.0. 
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-  The  aperture  has  an  x-dimensicn  of  0.15  cm  and  a  y-dimension  of 
1.1  cm . 

-  Both  patches  have  dimensions  of  2.5  x  2.5  cm  and  are  centered 
over  the  aperture. 

The  MoM  segmentation  results  in  45  x-directed  and  40  y— directed 
current  cells  on  each  patch.  On  each  patch,  a  -  0.25  cm  and  b  -  0.5  cm. 
The  MoM  impedance  matrix  is  170  x  170  and  each  impedance  sub— matrix  is 
85  x  85. 

5.1.2  Antenna  *r2 .  The  antenna  referred  to  as  antenna  #2  is 
identical  to  antenna  #1,  except  the  patch  of  interface  3b  is  2.65  x  2.65 
cm.  Measurements  of  this  design  shows  resonant  frequencies  for  this 
antenna  at  3.54  GHz  and  4.15  GHz  [30:51], 

The  MoM  impedance  matrix  and  sub-matrices  for  antenna  #2  have  the 
same  dimensions.  The  cell  dimensions  for  the  patch  on  interface  2b  are 
the  same  as  for  antenna  #1,  but  on  the  patch  of  interface  3b ,  a  *  0.265 
cm  and  b  -  0.53  cm. 


5.2  Surface  Current  Data 

Figures  5-1  through  5-10  are  plots  of  the  surface  current  data  as 
calculated  by  the  MoM  software.  The  plots  were  made  using  data  obtained 
from  the  files  JREAL.DAT  and  JIMAG.DAT  and  are  presented  in  pairs,  with 
the  real  component  of  the  surface  current  shown  in  the  left  plot  and  the 
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imaginary  in  the  right  plot.  As  discussed  in  Section  4.6.2,  subroutine 
JDISTR  produces  data  at  sample  points  located  in  the  center  of  the 
current  cells  along  the  row,  and  the  location  of  the  row  is  given  above 
each  plot. 

5.2.1  Surface  Currents  on  Antenna  #2.  Figure  5-1  shows  the  x- 
directed  currents  in  the  top  three  of  the  five  rows  of  x-directed 
current  cells  on  the  patch  of  interface  2b  for  antenna  #1.  Comparison 
of  the  magnitudes  of  the  real  and  imaginary  components  of  the  currents 
shows  the  imaginary  component  of  the  currents  are  clearly  dominant.  The 
reason  for  this  dominance  is  that  the  antenna  is  operating  near  reso¬ 
nance  at  3.7  GHz  and  at  resonance  the  surface  currents  are  always 
dominated  by  the  imaginary  component  [18:321].  Also,  both  components 
are  symmetric  along  the  x-  and  y-axes  w- ' '  aspect  to  the  cells  center. 

Figure  5-1  also  shows  how  the  real  components  of  the  currents  are 
fairly  constant  as  the  cell's  is  traversed  along  the  y-axis .  On  the 
other  hand,  the  imaginary  components  vary  appreciably  along  the  y-axis, 
peaking  in  magnitude  at  the  cell's  center  (see  plot  (c)  of  Figure  5-1). 
This  peak  in  the  imaginary  component  of  the  currents  means  the  strongest 
coupling  occurs  at  the  patch's  center.  These  same  features  are  evident 
in  the  x-directed  currents  on  the  patch  of  interface  3b  (see  Figure 
5-4). 

Figures  5-2  and  5-3  show  the  y-directed  currents  on  the  patch  of 
interface  2b.  Figure  5-2  shows  currents  in  the  first  three  rows  of 
y-directed  current  cells  on  the  patch,  and  Figure  5-3  shows  the  currents 
in  rows  four,  five,  and  six.  While  the  magnitudes  of  the  y-directed 
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Figure  5-2  Plots  of  y-directed  Surface  Currents  for  the  First  Three 
Rows  of  y-cells  on  the  Patch  of  Interface  2b,  Antenna  #1;  (data  from 
four  node  quadrature) 
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Figure  5-3  Plots  of  y-directed  Surface  Currents  for  Rows  Four,  Five, 
and  Six  of  y-cells  on  the  Patch  of  Interface  2b,  Antenna  #1  (data  from 
four  node  quadrature) 
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currents  are  much  smaller  than  those  of  the  x-directed  currents ,  the 
dominance  of  the  imaginary  component  of  the  currents  is  evident.  In 
contrast  to  the  x-directed  currents,  the  y-currents  appear  to  show  no 
symmetry  at  all.  However,  close  comparison  of  plots  (b)  and  (c)  in 
Figure  5-2  shows  the  surface  currents  are  symmetric  along  the  cell’s  x- 
axis  with  respect  to  the  center  of  the  cell,  but  as  the  center  of  the 
cell  is  crossed,  the  currents  experience  a  180®  phase  shift.  This 
feature  is  also  evident  in  all  data  obtained  for  the  y-directed  cur¬ 
rents.  (see  Figure  5-6,  plots  (b)  and  (c)  and  Figure  5-10,  plots  (a)  and 
(b)). 

Figures  5-4  through  5-6  show  the  surface  currents  on  the  patch  of 
interface  3b.  These  currents  have  the  same  general  features  as  those 
explained  for  the  patch  of  interface  2b.  As  compared  to  the  currents  of 

interface  2b,  the  currents  of  interface  3b  are  about  one  order  of 

magnitude  smaller  than  those  of  interface  2b.  Also,  a  comparison  of 
plots  (a)  and  (b)  of  Figures  5-1  and  5-4  shows  that  there  is  a  phase 

shift  of  180°  between  the  imaginary  components  of  the  currents  on  the 

two  patches.  As  was  the  case  with  the  patch  of  interface  2b,  the 
imaginary  component  of  the  currents  peaks  at  the  cell's  center. 

Figures  5-7  and  5-8  are  plots  of  the  same  currents  as  shown  in 
Figures  5-1  and  5-4,  respectively,  except  data  plotted  in  Figures  5-7 
and  5-8  were  obtained  using  six  node  quadrature.  While  the  distri¬ 
butions  of  the  real  current  components  have  the  same  general  shapes  for 
four  and  six  node  quadrature,  the  distributions  of  the  imaginary 
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Figure  5-4  Plots  of  x-directed  Currents  for  the  First  Three  Rows  of 
x-cells  on  the  Patch  of  Interface  3b,  Antenna  #1  (data  from  four  node 
quadrature) 
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Figure  5-5  Plots  of  y-directed  Currents  for  the  First  Three  Rows  of 
y-cells  on  the  Patch  of  Interface  3b,  Antenna  #1  (data  from  four  node 
quadrature) 
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Figure  5-6  Plots  of  y-directed  currents  for  the  First  Three  Rows  of 
y-cells  on  the  Patch  of  Interface  3b,  Antenna  #1  (data  from  four  node 
quadrature) 
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Figure  5-7  Plots  of  x-directed  Currents  for  the  First  Three  Rows  of 
x-cells  on  the  Patch  of  Interface  2b,  Antenna  #1  (data  from  six  node 
quadrature) 
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components  vary  significantly.  The  disagreement  between  results  is  best 
shown  by  a  comparison  of  plot  (c)  in  Figures  5-7  and  5-8.  This  disagree¬ 
ment  could  be  caused  by  an  inaccuracy  of  the  Green's  function  approxima¬ 
tions.  Tests  of  the  approximations  showed  the  approximations  were 
accurate  for  R  >  0.0008.  For  six  node  quadrature  with  source  and 
observer  located  in  the  same  cell,  the  values  of  R  can  go  well  below 
this  limit,  so  there  are  some  errors  in  the  Green's  function  for  this 
case.  Since  the  self-terms  of  the  impedance  matrix  are  dominant,  the 
error  in  the  Green's  functions  has  a  pronounced  impact  on  the  results. 

5.2.2  Surface  Currents  on  Antenna  # 2 .  Figure  5-9  shows  the  x- 
directed  currents  for  the  first  three  rows  of  x-cells  on  the  patch  of 
interface  3b  for  '  *  nna  #2.  Comparisons  with  Figure  5-4  show  the 
current  distr _h  .ion  are  generally  the  same,  but  the  current  magnitudes 
on  antenna  #2  are  larger  than  those  of  antenna  #1.  This  increase  in 
current  magnitude  is  most  likely  due  to  the  fact  that  the  larger  patch 
of  antenna  #2  results  in  stronger  coupling  with  the  sources  of  the 
antenna.  Comparisons  of  tabulated  current  distributions  data  show  the 
surface  current  magnitudes  on  the  patch  of  interface  2b  in  antenna  #2 
are  about  10%  smaller  than  the  same  currents  in  antenna  #1.  This 
finding  also  supports  the  idea  that  the  larger  patch  on  interface  3b 
results  in  stronger  electromagnetic  coupling  with  the  sources  in  the 
antenna . 

Figure  5-10  shows  the  y-directed  currents  in  rows  five  and  six  on 
the  patch  of  interface  3b  in  antenna  #2 .  Compared  to  the  same  y- 
directed  currents  on  antenna  #1,  these  currents  are  slightly  greater  in 
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Figure  5-9  Plots  of  the  x-directed  Currents  for  the  First  Three  Rows  of 


x-cells  on  the  Patch  of  Interface  3b,  Antenna  #2  (data  from  four  node 


quadrature) 
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magnitude,  but  essentially  identical.  Comparing  plots  (a)  and  (b)  of 
Figure  5-10  also  shows  the  y-directed  currents  of  antenna  #2  have  the 
same  symmetry  and  phase  shift  characteristics  as  antenna  #1. 
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5.3  Surface  Current  Expansion  Coefficients 


A  listing  of  the  surface  current/charge  density  expansion  coeffi¬ 
cients  obtained  with  four  node  Gaussian  numerical  integration  for  the 
patch  on  interface  2b  of  antenna  #1  is  given  in  Table  5-1  and  for 
interface  3b  is  given  in  Table  5-2.  Several  observations  concerning  the 
coefficients  can  be  made, which  are  valid  for  both  the  four  node  and  six 
node  solutions. 

As  illustrated  by  the  surface  current  plots  of  section  5.2,  the 
expansion  coefficients  for  the  x-directed  surface  currents  are  symmetric 
about  the  center  of  the  cell.  An  examination  of  x-directed  coefficients 
1  through  9  on  either  patch  shows  the  symmetry  in  x,  and  comparing 
coefficients  1  through  9  with  coefficients  37  through  45  shows  the 
symmetry  in  y. 

As  expected,  the  imaginary  component  of  the  coefficients  are 
dominant  over  the  real  component  on  both  patches.  On  interface  2b  the 
magnitude  of  the  imaginary  components  are  100  to  lO*1  times  larger  than 
those  of  the  real  components.  On  the  patch  of  interface  3b,  the 
imaginary  components  are  10  to  100  times  larger  than  the  real  compo¬ 
nents  .  — 
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Table  5-1  Surface  Current  Coefficients  for  the  Patch  of  Interface  2b 


1 — 

— — —  Hill 

ILLS  1 

INDEX  COEFFICIENT 

INDEX  COEFFICIENT 

1 

(0. 144930D-05 ,0. 313567D-02) 

24 

(-. 704937D-04 , - . 152963D-01) 

2 

(-. 588259D-04 ,-. 111001D-01) 

25 

(-. 490389D-04 , 0 . 550930D-02) 

3 

(-. 277699D-04 ,0. 389851D-01) 

26 

(-. 844589D-04 622061D-02) 

4 

(-.447774D-04 226170D-01) 

27 

(-. 1978800-04, 0.243996D-01) 

5 

(-. 390199D-04, 0.234268D-02) 

28 

(- . 180006D-04 , 0 . 733650D-02) 

6 

(-.447774D-04 226170D-01) 

29 

(-. 788197D-04 , 0 . 165448D-01) 

7 

(-. 277699D-04 ,0. 389851D-01) 

30 

(-.448388D-04 886093D-02) 

8 

(-.588259D-04.-.111001D-01) 

31 

(-.655276D-04,-. 5832160-02) 

9 

(0 . 144930D-05 , 0 . 313567D-02) 

32 

(-.5880500-04 , 0 . 156866D+00) 

10 

(- . 180006D-04 , 0 . 733650D-02 ) 

33 

(-. 655276D-04 583216D-02 ) 

11 

(-. 788197D-04 ,0. 165448D-01) 

34 

(- .448388D-04 ,-. 886093D-02 ) 

12 

-.448388D— 04 886093D-02) 

35 

(-. 788197D-04 ,0. 165448D-01) 

13 

<-. 655276D-04 , 583216D-02) 

36 

(-. 180006D-04 ,0. 733650D-02) 

14 

(-. 588050D-04.0. 156866D+00) 

37 

(0 . 144930D-05 , 0 . 313567D-02) 

15 

(-. 655276D-04 , - . 583216D-02) 

38 

(— . 588259D-04 ,-. 111001D-01) 

16 

(- . 448388D-04 . 88609 3D-0 2) 

39 

(— . 277699D-04 ,0 . 389851D-01) 

17 

(- . 788197D-04 , 0 . 165448D-01) 

40 

(-.447774D-04 ,-. 226170D-01) 

18 

(-. 180006D-04 , 0 . 733650D-02) 

41 

(-.390199D-04, 0.2342680-02) 

19 

(-.197880D-04 ,0.2439960-01) 

42 

(— .447774D-04 ,-. 226170D-01) 

20 

(-. 844589D-04 622061D-02) 

43 

(- . 277699D-04 , 0 . 389851D-01 ) 

21 

(- . 490389D-04 , 0 . 550930D-02) 

44 

(-. 588259D-04 ,-. 111001D-01) 

22 

(-. 704937D-04,-. 152963D-01) 

45 

(0 . 144930D-05 , 0 . 313567D-02 ) 

23 

(- . 635182D-04 , 0 . 410874D+00) 

Y-DIRECTED 

CURRENT 

CELLS 

1 

(0 . 309850D-05 , — . 126455D-04) 

21 

(0. 3 14 5 2 8D— 06 ,-. 271924D-04) 

2 

(-.434617D-06 600874D— 04) 

22 

(-.484581D-06 , 0 . 302674D-03) 

3 

(-. 276181D-05 872770D-04) 

23 

(0.650785D-06, 0.3772490-03) 

4 

(0 . 39881 8D-05 , 0 . 101 3 8 ID— 04) 

24 

(-. 155509D-06 ,-. 774813D-04) 

5 

(-. 292500D-06 . 674216D-04) 

25 

(0 . 152025D-06 ,- . 136765D-04) 

6 

(-. 186479D-05 . 365742D-04) 

26 

(-. 630530D-07 , 0 . 314760D-03 ) 

7 

(0 . 989116D— 06 . 577124D-04) 

27 

(0 . 717857D-06 , 0 . 389302D-03) 

8 

(-. 143753D-05 688197D-04) 

28 

(-. 174532D-06 ,-. 195733D-04) 

9 

(- . 893336D-07 306283D-04) 

29 

(0 . 119948D-06 , 0 . 233631D-04) 

10 

(-. 389707D-06 124038D-03) 

30 

(0 . 346488D-06 , 0 . 174468D-03) 

11 

(-. 346438D-06 174468D-03) 

31 

(0 . 389658D-06 , 0 . 124038D-03) 

12 

(-. 119989D-06 233631D-04) 

32 

(0. 893715D-07 , 0. 306283D-04) 

13 

(0. 174509D-06 , 0 . 195733D-04) 

33 

(0. 143755D-05 , 0 . 688197D-04) 

14 

(-. 717827D-06 389302D-03) 

34 

(- . 989152D-06 , 0 . 577124D-04) 

15 

(0 . 630119D-07 ,-.3147600-03) 

35 

(0 . 1864830-05 , 0 . 365742D-04) 

16 

(- . 151986D-06 , 0 . 136765D-04) 

36 

(0 . 292481D-06 , 0 . 674216D-04) 

17 

(0 . 155539D-06 ,0 . 774813D-04) 

37 

(-. 398819D-05 ,-. 101381D-04) 

18 

(-.650930D-06,-. 3772490-03) 

38 

(0 . 276182D-05 , 0 . 872770D-04) 

19 

(0. 4847370-06, -. 302674D-03) 

39 

(0.4346100-06,0. 600874D-04) 

20 

(-. 314572D-06 , 0 . 271924D-04) 

40 

(-. 309850D-05 , 0 . 126455D-04) 
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Table  5-2  Surface  Current  Coefficients  for  the  Patch  of  Interface  3b 


r— 

— 1  II  M  1  III  ■ 

1  INDEX  COEFFICIENT 

INDEX  COEFFICIENT 

1 

(0 . 502501D-03 ,-. 920439D-03) 

24 

(0. 292374D-03 , 0 . 157864D-01) 

2 

(0 . 202992D-03 ,0. 983445D-03) 

25 

(0.3301070-03 , 0. 853381D-02) 

3 

(0 . 306094D-03 ,-.1601720-01) 

26 

(0. 211369D-03 , 0 . 170281D-02) 

4 

(0. 2712130-03 .0.219468D-02) 

27 

(0. 5667110-03,-. 8651550-02) 

5 

(0 . 290519D-03 .  563160I>-02) 

28 

(0. 5312740-03 ,-. 3449920-02) 

6 

(0. 271213D-03 , 0. 219468D-02) 

29 

(0. 194061D-03 ,- . 675561D-02) 

7 

<0. 306094D-03 160172D-01) 

30 

(0.3067340-03,0.5876380-02) 

8 

(0 . 202992D-03 , 0 . 983445D-03 ) 

31 

(0 . 2720540-03 ,- . 180131D-02) 

9 

(0. 502501D-03 920439D-03) 

32 

(0.2887590-03, -.2364120-01) 

10 

(0. 531274D-03 344992D-02) 

33 

(0.2720540-03, -.1801310-02) 

11 

(0.1940610-03, -.67 55610-02) 

34 

(0 . 3067340-03 , 0 . 587638D-02) 

12 

(0 . 306734D-03 ,0 . 587638D-02) 

35 

(0 . 1940610-03 ,-. 675561D-02) 

13 

(0 . 272054D-03 . 180131D-02 ) 

36 

(0 . 5312740-03 ,-. 344992D-02) 

14 

(0 . 288759D-03 236412D-01) 

37 

(0 . 502501D-03 ,-. 920439D-03) 

15 

(0 . 272054D-03 180131D-02 ) 

38 

(0. 2029920-03, 0 . 983445D-03) 

16 

(0 . 306734D-03 , 0 . 587638D-02) 

39 

(0 . 306094D-03 ,- . 160172D-01) 

17 

(0 . 194061D-03 675561D-02) 

40 

(0 . 271213D-03 , 0 . 219468D-02) 

18 

(0 . 531274D— 03 344992D-02 ) 

41 

(0 . 290519D-03 ,- . 563160D-02 ) 

19 

(0 . 566711D-03 86515 5D-02) 

42 

(0 . 271213D— 03 , 0 . 219468D-02 ) 

20 

(0 . 211369D-03 , 0 . 170281D-02) 

43 

(0. 306094D-03 160172D-01) 

21 

(0 . 330107D-03 , 0 . 853381D-02) 

44 

(0 . 202992D-03 , 0 . 983445D-03) 

22 

(0 . 292374D— 03 , 0 . 157864D-01) 

45 

(0 . 502501D-03 . 920439D-03) 

23 

(0 . 310082D-03 , 0 . 2 28 84 7D— 01) 

Y-DIRECTED 

CURRENT  CELLS 

1 

(0 . 288394D— 04 , 0 .480673D-03) 

21 

(-. 195355D-05 177182D-03) 

2 

(-.484272D-05 ,0 .450811D-04) 

22 

(0. 521850D-07 ,-. 201216D-02) 

3 

(-. 116127D-04 158249D-03) 

23 

(0 . 141238D-05 , 216219D-02) 

4 

(0 . 315246D-04 , 0 . 555294D-03) 

24 

(- . 248727D-05 ,— . 360807D— 03) 

5 

(0. 157141D-04.0. 599795D-03) 

25 

(-. 581689D-05 ,— . 101351D-02) 

6 

(-. 805477D-05 232144D-05) 

26 

(0. 106620D-05,-. 151859D-02) 

7 

125982D-05 45606 5D-04) 

27 

(0. 312257D-05 ,-. 107837D-02) 

8 

(0. 128715D-04.0.646055D-03) 

28 

(-. 661839D-05 106295D-02) 

9 

(0 . 107619D-04 , 0 . 88662 ID-03) 

29 

(-.1025650-04,-. 9792070-03) 

10 

(—. 4005 50D— 05,0. 173890D-03) 

30 

(0 . 274035D-05 . 178920D-03) 

11 

(-. 274037D-05 , 0 . 178920D-03) 

31 

(0 . 400551D-05 173890D-03) 

12 

(0. 102565D-04 , 0 . 979207D-03) 

32 

(-.1076190-04,-. 886621D-03) 

13 

(0 . 661841D-05 , 0 . 106295D-02) 

33 

(-.128715D-04, -.6460550-03) 

14 

(-. 312259D-05 ,0 . 107837D-02 ) 

34 

(0 . 125983D-05 , 0 . 456065D-04) 

15 

(-. 106618D-05 , 0 . 151859D-02) 

35 

(0 . 805476D-05 , 0 . 232144D-05) 

16 

(0 . 581689D-05 ,0 . 101351D-02) 

36 

(-.1571400-04,-. 5997950-03) 

17 

(0 . 248729D-05 , 0 . 360807D-03) 

37 

(-. 315246D-04 ,-. 5552940-03) 

18 

(-. 141244D-05 ,0. 216219D-02) 

38 

(0. 116127D-04 ,0. 158249D-03) 

19 

(-.5211660-07,0. 201216D-02 ) 

39 

(0.484272D-05,-. 4508110-04) 

20 

(0. 1953510-05, 0.177182D-03) 

40 

(-.288394D-04,-. 4806730-03) 
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5.4  MoM  Impedance  Matrix  Properties 


Due  to  space  limitations,  the  complete  HoM  impedance  matrix  cannot 
be  included  as  part  of  this  document,  however  the  matrix  elements  were 
examined  and  this  section  reports  some  of  the  observed  properties  of  the 
matrix. 

As  discussed  in  section  3.3.4,  the  MoM  impedance  matrix  is 
expected  to  be  diagonally  dominant,  and  the  matrix  produced  by  the  MoM 
software  did  show  the  dominance  of  the  diagonal  elements.  Additionally, 
as  is  the  case  at  resonance  [20:428],  the  imaginary  components  of  the 
elements  were  dominant  over  the  real  components.  Also,  the  elements 
describing  x-directed  fields  created  by  x-directed  sources  were  at  least 
one  order  of  magnitude  larger  than  the  remaining  elements  -  i.e  those 
describing  y-directed  fields  created  by  x-directed  sources.  This 
dominance  manifests  itself  in  the  relatively  large  values  of  the 
coefficients  for  the  x-directed  cells  as  compared  with  the  coefficients 
for  the  y-directed  cells. 

Each  impedance  sub-matrix  is  composed  of  four  smaller  impedance 
sub-sub-matrices  arranged  as  shown 


(5-1) 
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where  the  first  letter  of  the  superscript  is  the  observer  orientation 
and  the  second  is  the  source  orientation.  The  number  of  rows  in  the 
sub-sub-matrix  are  determined  by  the  number  cells  with  the  observer 
orientation  and  the  number  of  columns  is  determined  by  the  number  of 
cells  with  the  source  orientation.  For  instance,  consider  antenna  #1 
with  45  x-current  cells  and  40  y-current  cells  on  the  patch  of  interface 
2b.  Therefore,  the  sub-sub-matrix  with  superscript  xy  will  have  the 
dimension  45  x  40. 

A  study  of  the  first  45  x  45  block  of  the  impedance  matrix, 
representing  the  elements  describing  x-directed  fields  created  by  x- 
directed  sources  on  interface  2b,  also  shows  several  useful  symmetry 
characteristics.  First,  the  entire  45  x  45  block  consists  of  five, 
basic  9x9  sub-sub-matrices  rearranged  in  a  simple  pattern.  There  are 
five  basic  sub-sub-matrices  because  there  are  five  rows  of  x-directed 
cells  on  the  patch.  The  dimension  of  each  basic  sub-sub-matrix  is  equal 
to  the  number  of  x-directed  current  cells  on  the  patch.  If  the  elements 
in  rows  1  through  9  and  columns  1  through  9  are  designated  as  sub-sub¬ 
matrix  zn,  and  the  elements  in  rows  1  through  9  and  columns  10  through 
18  are  designated  as  z12,  and  so  forth  so  that  the  last  sub-sub-matrix 
comprised  by  elements  in  rows  1  through  9  and  columns  37  through  45  is 
designated  z15,  the  remaining  twenty  9x9  blocks  of  the  45  x  45  block 
can  be  filled  according  to  the  formula 


Z»,n  "  Zl, |m-n | 


(5-2) 
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where 


m  -  the  row  Index  for  the  sub-sub-matrix 
n  -  the  column  for  the  sub-sub-matrix 

This  pattern  is  the  same  as  that  of  the  toeplitz  MoM  impedance 
matrix  for  a  wire-scatterer ,  except  9x9  sub-sub-matrices  are  rear¬ 
ranged  instead  of  individual  matrix  elements.  Each  9x9  sub-sub-matrix 
is  toeplitz,  and  the  element  arrangement  pattern  is  given  by  equation 
(5-2)  with  the  z's  being  the  individual  matrix  elements  in  the  sub-sub- 
matrix . 

Since  the  symmetry  of  the  matrix  elements  is  a  due  to  the  fact 
that  the  Green's  functions  are  dependent  only  upon  the  argument  Jci?,  it 
is  sensible  to  believe  the  three  remaining  sub-sub-matrices  of  Zb22  are 
symmetric  and  toeplitz.  This  symmetry  is  important  because  it  can  be 
used  to  reduce  the  computer  processor  time  needed  to  fill  the  MoM 
impedance  matrix.  For  instance,  with  the  symmetry  explained  for  the 
first  45  x  45  matrix,  the  block  could  be  filled  by  first  calculating  and 
then  rearranging  the  45  elements  of  the  first  row.  Based  on  this 
observation,  it  is  sensible  to  expect  that  the  complete  MoM  matrix  can 
be  filled  by  calculating  only  four  rows  from  the  matrix.  Calculating 
the  first  row  would  provide  the  elements  needed  to  fill  the  sub-sub¬ 
matrices  for  x-directed  fields  on  interface  2b  that  ere  created  by 
sources  on  both  patches.  Calculating  the  row  with  an  index  of  (N2jc  +  1) 
would  provide  the  elements  needed  to  fill  the  sub-sub-matrices  describ¬ 
ing  y-directed  fields  on  interface  2b  created  by  sources  on  both 
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patches.  Data  for  the  row  with  index  (N2  +  1)  would  provide  the  ele¬ 
ments  needed  to  fill  the  sub-sub-matrices  describing  x-directed  fields 
on  interface  3b  created  by  sources  on  both  patches.  Finally,  elements 
with  a  row  index  equal  to  (N2  +  N3x  +1)  can  be  used  to  fill  the  sub¬ 
sub-matrices  describing  the  y-directed  fields  on  interface  3b  created  by 
the  sources  on  both  patches . 

Modifying  the  MoM  solution  software  to  exploit  the  toeplitz  nature 
and  symmetry  of  the  impedance  matrix  would  be  very  easy  to  do  and  would 
save  considerable  computer  processor  time.  While  the  exact  savings  is 
hard  to  guess  without  performing  the  modifications  and  running  the 
software,  estimates  can  be  made.  For  instance,  consider  the  solution 
for  antenna  #1  which  has  a  170  x  170  impedance  matrix  or  28,900  imped¬ 
ance  elements.  Calculating  only  four  rows  of  the  impedance  matrix  would 
require  finding  only  680  elements.  Since  each  element  calculation 
involves  the  execution  of  a  three-layer,  nested  do-loop,  -  one  loop  for 
the  line  integration  and  two  for  the  surface  integration  -  it  is 
reasonable  to  believe  most  of  the  processor  time  for  the  program  is 
dedicated  to  filling  the  impedance  matrix.  Since  reducing  the  number  of 
element  calculations  from  28,900  to  680  is  a  97%  reduction,  a  similar 
reduction  in  the  processor  time  should  be  obtained. 
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VI.  Recommendations  and  Conclusions 


This  chapter  contains  conclusions  about  the  accuracy  of  the  MoM 
solution  developed  in  this  effort  and  recommendations  for  improving  the 
MoM  software's  flexibility  and  capabilities. 


6.1  Conclusions 

Based  on  data  presented  in  Chapter  V,  the  data  obtained  from  the 
MoM  solution  software  is  correct,  as  best  as  can  be  told  at  this  point. 
As  discussed  in  Section  5.2.1,  the  data  from  the  MoM  software  predicts 
the  antennas  are  near  resonance  at  3.7  GHz,  and  this  result  matches 
measurements.  The  solution  predicts  the  strongest  currents  occur  in  the 
center  of  the  patches  which  is  expected  since  both  patches  are  aligned 
directly  over  the  aperture  of  the  antenna.  Also,  The  MoM  solution 
predicts  a  dominance  of  the  x-directed  surface  currents  over  the  y- 
directed  currents,  and  this  result  is  expected  since  the  magnetic 
surface  currents  of  the  aperture  are  in  the  y-direction  only. 
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Unfortunately,  at  the  present  time  the  ability  to  derive  impedance 
information  and  field  patterns  using  data  from  the  MoM  software  has  not 
been  developed.  Since  this  data  is  the  standard  by  which  solution 
accuracy  is  judged,  the  final  measure  of  the  accuracy  of  the  data  from 
MoM  software  cannot  be  made. 


6.2  Expanding  the  MoM  Software 

Microstrip  antenna  analysis  software  is  most  commonly  used  for 
predicting  the  antenna  impedance  and  radiation  characteristics.  At 
present,  the  MoM  solution  software  is  not  capable  of  providing  data  for 
either  predictions.  This  section  outlines  the  modifications  and 
additions  that  could  be  made  to  use  the  MoM  software  for  impedance  and 
field  pattern  predictions. 

6.1.1  Incorporating  a  Feed  System  Model.  Before  any  antenna 
impedance  data  can  be  calculated,  the  MoM  software  must  have  access  to 
data  describing  the  feed  system  operation.  With  this  data,  the  complete 
antenna  model  of  Sections  3.2  and  3.3  can  be  completed,  thus  providing  a 
means  to  completely  calculate  the  scattered  and  excitation  fields  on  the 
feedline.  Once  the  excitation  and  scattered  fields  are  known,  the  input 
impedance  can  be  calculated  from  the  relationship  [20:428] 
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2in  -  -^Joh(£<l>  *£'•>)•  idz 


(6-1) 


where 

I  -  the  Input  current 

h  -  the  distance  between  the  ground  plane  and  the  feedline 

£(,)  -  the  scattered  fields  on  the  feedline 
.  the  excitation  fields  on  the  feedline 

Once  input  impedance  can  be  calculated,  VSWR  and  other  performance 
data  can  be  obtained.  Additionally,  the  MoM  impedance  matrix  can  be 
analyzed  to  obtain  data  on  resonant  frequencies  of  the  antenna. 

It  should  also  be  noted  that  the  completion  of  the  MoM  solution 
for  the  full  model  would  require  only  minor  modif ications  to  the  MoM 
software,  however  the  computer  memory  and  processor  time  requirements 
would  increase  appreciably. 

6.1.2  Calculations  for  Bandwidth  Characteristics .  Calculating 
the  antenna  characteristics  over  a  bandwidth  requires  the  ability  to 
evaluate  the  Green's  functions  for  various  frequencies.  There  are 
several  ways  this  capability  might  be  incorporated  into  the  MoM  soft¬ 
ware  . 

The  simplest  option  would  require  creating  Green's  function 
polynomial  coefficient  files  for  selected  sampling  frequencies  over  the 
bandwidth  of  interest.  As  the  frequency  of  the  excitation  is  varied, 
the  MoM  software  would  read  the  polynomial  coefficients  for  the  present 
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frequency,  and  then  complete  the  MoM  solution  for  the  frequency. 

However,  incorporating  this  option  would  require  all  Green’s  functions 
polynomials  to  have  the  same  general  form  for  all  frequencies  of 
interest  which  is  probably  not  the  case. 

Other  ways  to  enable  bandwidth  studies  are  using  Green's  functions 
polynomials  that  are  a  function  of  both  R  and  frequency,  incorporating 
the  Green's  functions  evaluation  code  into  the  MoM  software,  or  imple¬ 
menting  an  interpolation  procedure  using  tabulated  data. 

These  three  options  would  all  result  in  increased  processor  time 
for  the  MoM  solution.  Using  frequency  dependent  polynomial  coefficient 
data  is  the  most  efficient  means  in  terms  of  computer  memory.  Using  an 
interpolation  scheme  would  significantly  increase  memory  requirements 
due  to  storage  for  the  tabulated  Green's  function  data.  Incorporating 
the  code  to  calculate  the  Green's  functions  would  produce  the  most 
accurate  results,  but  the  accuracy  improvement  would  most  likely  be 
offset  by  an  enormous  increase  in  computation  time.  However,  the 
additional  computer  time  would  be  almost  identical  to  the  time  needed  to 
produce  the  data  needed  for  interpolation  tables  or  to  produce  polynomi¬ 
al  approximations. 

6.1.3  Modifications  for  Field  Pattern  Calculations .  Calculating 
the  field  patterns  of  the  antenna  could  easily  be  incorporated  into  the 
MoM  solution  software.  The  major  obstacle  to  this  capability  is 
obtaining  the  closed  form  expressions  for  the  gradient  of  the  Gq  Green's 
functions.  Without  these  expressions,  the  fields  created  by  the  surface 
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charge  distributions  cannot  be  evaluated  numerically,  and  therefore  the 
fields  on  the  antenna  interfaces  cannot  be  calculated. 

With  the  gradient  of  the  Gq  functions  known,  a  subroutine  to 
establish  the  observer  location  and  calculate  the  contributions  from  the 
aperture  and  patches  of  the  antenna  would  need  to  be  written.  The 
functions  performed  by  the  field  calculation  subroutine  would  be  exactly 
like  those  performed  by  subroutine  EVALSUB,  with  the  exception  that  the 
integration  along  the  observer  test  segment  would  not  be  performed. 

Once  the  fields  of  the  antenna  interface  can  be  calculated,  the  software 
to  complete  the  calculations  of  the  aperture  model  as  explained  in 
Section  3.7  would  need  to  be  written  and  incorporated  into  the  MoM 
solution  software.  With  these  additions,  the  field  patterns  of  the 
antenna  could  be  calculated. 


6.3  Improving  Computation  Efficiency 

As  mentioned  in  Section  5.4,  the  3-plane  MoM  impedance  matrix  is 
toeplitz  and  symmetrical.  While  the  symmetry  is  complex,  the  reduction 
in  the  number  of  impedance  elements  obtained  by  taking  advantage  of  the 
symmetry  and  toeplitz  nature  of  the  matrix  would  make  the  study  of  the 
symmetry  worthwhile. 
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Appendix  A.  Green's  Functions  Somaerfeld  Integrals  and 
Polynomial  Approximations 


This  appendix  lists  the  Sommerfeld  integrals,  as  derived  by  Nazar, 
that  define  the  Green's  functions  for  the  aperture-fed,  stacked  patch 
antenna.  For  the  case  when  the  observer  and  source  plane  coincide,  the 
asymptotic  representation  of  the  integral  for  R  -*  0  is  also  given. 

These  integrals  are  taken  directly  from  Nazar's  thesis  [23:112-122]. 
Additionally,  the  coefficients  for  the  Green's  function  polynomial 
models,  as  derived  by  Erwert  and  used  in  the  MoM  solution,  are  tabulat¬ 
ed. 

A.l  Definitions  of  Constants  and  Variables  and  Polynomial  Model  Forms 

In  all  of  the  Sommerfeld  integrals  the  path  of  integration  C  is 
along  the  real  axis  in  the  complex  A-plane,  except  for  deformations 
around  any  integrand  singularities.  The  deformations  around  the  poles 
must  be  made  in  the  half-plane  such  that  Im(A)  >  0  [21:50].  Following 
are  the  definitions  of  the  constants  used  in  the  Green's  function 
Sommerfeld  integrals. 
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Mib’  M 2b.  P3b  “  the  permeability  of  the  dielectric  between  inter¬ 
face  lb  and  2b,  interface  2b  and  3b,  and  above  interface  3b, 
respectively;  p3b  is  actually  /i0 

Mb23  “  the  ratio  Mzt/Mab-  all  other  subscripted  p  with  two  numbers 
in  the  subscript  are  defined  in  the  same  way 

blb,  hjt  -  the  distance  between  interfaces  lb  and  2b,  and  interfac¬ 
es  2b  and  3b,  respectively 

H0(2)(AR)  -  the  Hankel  function  of  the  second  kind  of  order  0  with 
argument  AR 


The  polynomial  models  for  the  Green's  functions  consists  of 
separate  models  for  the  real  and  imaginary  components.  For  the  GA 
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Green's  functions  and  the  Ir,  I*,  and  1 1  components  of  the  GE  Green's 


functions,  the  real  component  is  represented  by  an  eighth  order  *yno- 
mial  for  0.001  <  R  <  0.01.  When  R  >  0.01,  the  real  component  is  modeled 
as  follows: 


c  c 

Imag{  f(R)  )  -  +  _i  +  c2  +  c3R  +  ckR2  +  csR3  + 

c6R*  *  c7R 5  +  caR6  *  ceR7  +  c10Re 


(A-2) 


where  the  function  f(R)  is  either  one  of  the  GA  functions  or  one  of  the 
components  of  the  GE  functions. 

The  imaginary  components  of  the  GA  functions  and  the  components  of 
the  Ge  functions  are  represented  by  a  sixth  order  polynomial  valid  for 
all  values  of  R. 

The  real  components  of  the  Gq  functions  are  modeled  by  a  sixth 
order  polynomial  valid  for  all  R.  The  imaginary  components  are  modeled 
by  an  eighth  order  polynomial  for  0.001  <  R  <  0.01,  and  by  a  polynomial 
as  given  in  equation  (A-2)  for  values  of  R  >  0.01. 

The  coefficients  for  each  polynomial  model  are  tabulated  along 
with  the  Sommerfeld  integral  defining  the  Green's  functions.  In  the 
cases  where  a  component  of  the  function  is  represented  by  two  different 
polynomials,  the  coefficients  listed  under  the  column  heading  "Real  #1," 
for  example,  would  indicate  the  coefficients  for  the  polynomial  model - 
when  R  «  [0.001,  0.001]  and  the  column  heading  "Real  #2"  would  indicate 
the  coefficients  for  the  polynomial  model  when  R  >  0.001. 
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A.  2  Green's  Functions  Describing  Fields  Generated  by  Sources  on 
Interface  lb. 


A. 2.1  Green's  Functions  for  Fields  on  Interface  2b. 


<?&?(*, D  --sin(2f)  ^Ibai(«)-Ib2iw] 
clliiR.O  -  -[-Ib2i<«) -coS2(f)Ibr21(tf)  +  cos<2f)Ib821(tf)] 
GtllUl,!)  -  [-l£2100  +sin2<f)l£21(*)  +cos(2f)Ib21w] 
ctll(R,n  -  -Sin(2f)  [-  1i£21(K)  +ll21(R) 


where 


«■<*>  ‘  J, 


h£2)(\R)  sech  (  blbula)  A: 


- - r2(A)  dA 

>*(• 

Jc  RD*(\) D.( X) 


c  D.°(A)D°(A) 

(2), 

V' 

“rd:(A)D°(A) 

(2), 


rt  r  #0 ^(AJ?)  sech(blbula)  A 

•  bZl  “  L - r - 

Jc  d*<a) 


£bi3u3bC°s^(u2b  (b2b  -  blb) ) 
+  ebl2u2bsinh(u2b  (^2b  "  ^ lb)  ) 


and 


r2(A)  -  (  1  -  Mbl3  £bl3  )  u 


2b 


(  -  —  Pbl3  £  b  13  )  u2b 

2 

*  (  ~Mbl3  ebl3  +  Mb23  £b23  )  u3b 


sinh^UjtCbjb-b^)] 


^2b  ^3b 


^b23  "  Mbl3£bl2 
~A‘bl2£bl3  +  £  b23 


u2bu3b 


sinh[2u2b  (b2b-blb) 


(A-3) 


(A-4) 


dA 


(A-5) 


Table  A-l  Coefficients  for  Irb2i  Polynomial  Model 


Coefficient  # 


Real  #1 


-5 . 1497773E+03 


2 . 0098000E+09 


1 . 1664054E+11 


-2.4628222E+14 


Real  #2 


Imaginary 


-2 . 3690010E+04  -2 . 7567323E-01  7 . 3622798E+02 


4 . 1969248E+01 


1 . 6413327E+02 


-1 . 7263071E+03  |  -9 .4660002E+05 


6 . 7013484E-04  -5 . 5624401E+06 


6 . 0581647E+16 


2 . 5393666E+06 


-1 . 0559782E+00 


-3 . 8054523E+09 


1 . 0860188E+11 


-1 . 3359528E+12 


7 . 9367587E+12 


-1 . 8720567E+13 


6 . 6599743E+08 


Table  A-2  Coefficients  for  Isb21  Polynomial  Model 


Table  A-3  Coefficients  for  lS>2i  Polynomial  Model 


Coefficient  # 

Real  #1 

Real  #2 

Imaginary 

0 

8 . 3730489E+05 

1 . 0488087E+00 

-1 . 2618585E+03 

1 

1 . 8548911E+07 

-9 . 7803738E+01 

-1 . 1645510E+01 

2 

-6 . 9596640E+11 

3 . 8139458E+03 

1 . 2324177E+06 

3 

4 . 9585244E+14 

-1 . 1916493E-03 

7 . 8552878E+06 

4 

-1 . 6509095E+17 

-4 . 5203333E+06 

-8 . 1658149E+08 

5 

3 . 0643030E+19 

1 . 8826109E+00 

9 . 7146319E+09 

6 

-3 . 2457653E+21 

6 . 2316114E+09 

-3.4550352E+10 

7 

1 . 8334913E+23 

-1.7672322E+11 

- 

8 

-4 . 2836938E+24 

2 . 1844777E+12 

- 

9 

- 

-1 . 3111428E+13 

- 

10 

- 

3 . 1297587E+13 

- 

-6 


A. 2. 2  Green's  Function  for  Fields  on  Interface  3b. 


cUliR.O  --sin(2f)  -£lbr31<iO -lZzi(R) 

Gf£(K,f)  -  -[-IbaiW  ~cos2(Dl£3i<K>  ♦cos(2f)l£3iW] 

G&y(J2,f)  -[-IbsiW  ♦sin2(f)Ibr31(JO+cos(2f)Ib‘31<io] 


(A— 6) 


GlUiR.O  -  -sin(20 


-j-^bsiw  *  -fb3i(£) 


rilbaiW^i 


where 


(R)  -  ( 

•».«>(»*)  X>  , 

J 

c  D*(A)d£(A) 

<*>-j 

'  4^ilrj(A)dx 

(A-7) 

c  i?D*(A)D*(A) 

(R)  -  I 

r  J»0(  ^  AlijJj  r  ,  ,  *  i  1  J, 

L— - c - —  ^bl3U3bsech(Ulb^lb)  ]dx 

J 

C  Ob (A) 

and 


r3(A)  -  (l-/ibl3«bl3)u!bsech(blbulb)cosh[u2b(b2b-blb)] 

(£b12-/ib23£bl3)ulbu2bCSC^(^lbulb) 

+ 

+  (  "^bU  £b!3  'f  £b23  )  u2bU3bSeC^  ^lbulb  ) 


(A-8) 
sinh[u2b  (b2b-blb )  ] 
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Table  A-4  Coefficients  for  Irb31  Polynomial  Model 


Coefficient  # 

Real  #2 

Real  #2 

Imaginary 

0 

-7 . 1631702E+04 

-3 . 9145822E-01 

1.4578382E+03 

1 

-1 . 6126649E+06 

6 . 0771449E+01 

3 . 2605113E+02 

2 

1.6138582E+10 

-2.4640424E+03 

-1 . 8752160E+06 

3 

-4 . 6090776E+12 

1 . 0734826E-03 

-1 . 1014941E+07 

4 

2 . 0514025E+14 

4.2372332E+06 

1 . 3194275E+09 

5 

1 . 1322934E+17 

-1 . 6926680E+00 

-1 . 6258111E+10 

6 

-2 . 2354787E+19 

-6 . 5280178E+09 

5 . 9202168E+10 

7 

1 . 6596165E+21 

1 . 8569636E+11 

- 

8 

-4 . 5322535E+22 

-2 . 2676061E+12 

- 

9 

- 

1 . 3351982E+13 

- 

10 

- 

-3 . 1213003E+13 

- 

Table  A-5  Coefficients  for  I*b31  Polynomial  Model 


Coefficient  # 

Real  #2 

Real  #2 

Imaginary 

0 

-3 . 5845677E+04 

-3 . 2220387E-01 

7 . 2947448E+02 

1 

-1 . 8303401E+05 

-1 . 1195931E+01 

-2 . 0927389E+02 

2 

3 . 3506986F ’ 09 

6 . 7963919E+02 

-4 . 3564658E+05 

3 

-4 . 0133292^+11 

-2 . 7539833E-04 

-3 . 3273908E+06 

4 

-1 . 6475295E+14 

-7 . 0817217E+04 

2.4440619E+08 

5 

5 . 8048203E+16 

4.4877070E-01 

-2 . 5345126E+09 

6 

-7 . 7186016E+18 

1 . 6953818E+08 

8 . 0937999E+09 

7 

4 . 8748677E+20 

-8 . 6162929E+09 

- 

8 

-1 . 2131994E+22 

1 . 5436197E+11 

- 

9 

- 

—1 . 1848456E+12 

- 

10 

- 

3 . 3469197E+12 

- 
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Table  A-6  Coefficients  for  1*1,31  Polynomial  Model 


Coefficient  # 

Real  #1 

Real  #2 

Imaginary 

0 

2 . 8114574E+05 

1 . 5538455E+00 

-2 . 4949310E+03 

1 

5 . 9078765E+06 

-1 . 2252199E+02 

-2 . 4719206E+01 

2 

-5 . 9241731E+10 

4. 5582187E+03 

2.4398514E+06 

3 

1 . 6899417E+13 

-1.4436500E-03 

1 . 5547730E+07 

4 

-7 . 7355313E+14 

-6.4603264E+06 

-1 . 6171597E+09 

5 

-4 . 0848676E+17 

2 . 2647814E+00 

1 . 9242115E+10 

6 

8 . 1118069E+19 

9 . 1346088E+09 

—6 . 8442981E+10 

7 

-6 . 0339466E+21 

-2.5611943E+11 

- 

8 

1 . 6495873E+23 

3 . 1131480E+12 

- 

9 

- 

-1 . 8348861E+13 

- 

10 

- 

4 . 3053899E+13 

- 
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A. 3  Green's  Functions  Describing  Fields  Generated  by  Sources  on 
Interface  2b. 

A. 3.1  Green's  Functions  for  Fields  on  Interface  2b. 


rbxx/p'l  - 
°K22  <K)  - 


r  « o2)< 

Jc  ~D? 


(X  R)X 


(A) 


u^coshfujtCbjb-bit,)] 
h  Mb23u3bs^-n^[u2b  (^2b  “  ^lb)  ] 


dX 


r*  byy 
JA22 


(JO 


^A22  (R) 


b  1  r  H<2)(XR)X 

q22  ‘  V  Dicx)Di(x)  q22 


dX 


(A— 9) 


where 


r 


q22 


£b23u2bu3bcos^2^  U2b  (^2b  ~^lb)  ) 


A*b23u2bu3b  + 


2  2 
(A»bl2u2b  +^bl3  eb23u3b) 

•UibtanhCbjbUjb) 


u2b  +  ^b 23  e b23  u3b  + 

+  < 

[^bl3  +^bl2  tb23]ulbu3btan^(^lbulb  ) 
+  ^bl3  eb23ulbu2btan^(  ^lbUlb  ) 


j  jsinh2[u2b  (b2b-blb)  ] 

sinh[2  u2b  (b2b-blb)  ] 


(A— 10) 


A.  3. 1.1  Asymptotic  Forms  for  R  -*  0 . 


Plb /Pe 

‘  T¥\Jo 


lL2(i?,A)dA  -fu 


PRAC  J0(X) 

Jo  T" 


(A-ll) 
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where 


T  -  1  *_  **b23 

^  1  +  Mbl2  +  J*b23  +  ^bl3 


(A-12) 


Gq22  (  R  ) 


where 


(A-13) 


1  *  t 


b23 


rq2 


(A-14 ) 


1  *  £b!2  +  eb23  +  ebl3 


and  the  integrals  I^CAR)  and  Ibq22(AR)  are  the  integrands  of  the 
original  Sommerfeld  integrals  with  H0(2)(AR)  replaced  by  J0(XR) . 
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Table  A-7  Coefficients  for  Polynomial  Model 


Coefficient  # 

Real  #1 

Real  #2 

Imaginary 

0 

2 . 8271357E-04 

1 . 1336337E-10 

-8 . 2697879E-08 

1 

-3 . 9846556E-01 

-7 . 7506745E-0 

3 . 0039534E-08 

2 

2 . 7285907E+02 

2 . 8775042E-07 

4.4481533E-05 

3 

-1.0962870E+05 

-7 . 9761966E-14 

3 . 7821366E-04 

4 

-4.2391776E+09 

-2 . 8113896E-04 

-2 . 3998044E-02 

5 

3 . 9943857E+11 

1 . 2917902E-10 

2 . 3993749E-01 

6 

-2 . 0833405E+13 

3.4524106E-01 

-7.4304513E-01 

7 

2 . 7278577E+07 

-9 . 6808158E+00 

- 

8 

4 . 6074276E+14 

1 . 2043431E+02 

- 

9 

- 

-7 . 3386638E+02 

- 

10 

1 . 7830531E+03 

- 

Table  A-8  Coefficients  for  Gbq22  Polynomial  Model 


Coefficient  # 

Real 

Imaginary  #1 

Imaginary  #2 

0 

2 . 4690865E-01 

-5 . 0295804E+02 

-3.0506451E-04 

1 

8 . 1564715E-02 

6 . 9944277E+05 

3 . 7573808E-02 

2 

-3 . 1310965E+02 

-4 . 7773874E+08 

-1 . 5014618E+00 

3 

-1 . 7501311E+03 

1 . 9210580E+11 

3 . 1525211E-07 

4 

2 . 1602490E+05 

-4.7861591E+13 

7 . 9630909E+02 

5 

-2 . 6605711E+06 

7.4454262E+15 

-5 . 0947594E-04 

6 

9 . 6760956E+06 

-7 . 0205947E+17 

-7 . 9604799E+05 

7 

- 

3 . 6635682E+19 

2 . 3725931E+07 

8 

- 

-8 . 1050716E+20 

-3 . 2364262E+08 

9 

- 

- 

2 . 1545771E+09 

10 

- 

- 

-5 . 6293092E+09 
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A. 3. 2  Green's  Functions  for  Fields  on  Interface  3b. 


where 

rq32  “  [u3b  +  Mbl3ulbtanh<blbulb)]u2bcosh[u2b(b2b-blb)] 

* 

+  Mbl3«bl2^b*Nb(Mb23-Mb13«bl2)]  .sinh[U2b(b2b-blb)] 

+  ^b!2ulbu3btan^(blbulb) 


(A-15) 


(A-16) 
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Table  A-9  Coefficients  for  GbA32  Polynomial  Model 


Coefficient  # 

Real  #1 

Real  #2 

Imaginary 

0 

4 . 3605698E-05 

1 . 7029449E-10 

-1 . 6312406E-07 

1 

1 . 6428591E-03 

-8 . 1399245E-09 

4 . 7418351E-08 

2 

-1 . 9408505E+01 

2 . 6622157E-07 

8 . 8848203E-05 

3 

1.0921513E+04 

-4 . 1338561E-14 

7 . 1341936E-04 

4 

-3 . 0476872E+06 

-2 . 9963275E-04 

-4 . 6862616E-02 

5 

4.9365870E+08 

6 . 3114878E-11 

4 . 6971581E-01 

6 

-4 . 6983776E+10 

3.4697986E-01 

— 1.4547107E+00 

7 

2 . 4369181E+12 

-9 . 2439916E+00 

- 

8 

-5 . 3137854E+13 

1 . 0942193E+02 

- 

9 

- 

-6 . 3847649E+02 

- 

10 

- 

1 . 4970486E+03 

- 

Table  A-10  Coefficients  for  Gbq32  Polynomial  Model 


Coefficient  # 

Real 

Imaginary  #1 

Imaginary  #2 

0 

4 . 8908078E-01 

-1 . 0079854E+02 

-4 . 6036133E-04 

1 

9 . 7940015E-02 

-3 . 9690570E+03 

5 . 4394079E-02 

2 

-6 . 1506941E+02 

4 . 6810350E+07 

-2 . 1168582E+00 

3 

-3 . 6456986E+03 

-2 . 6387382E+10 

2 . 8577567E-07 

4 

4 . 3101735E+05 

7 . 3656913E+12 

8.4397927E+02 

5 

-5 . 2961144E+06 

-1 . 1930563E+15 

-4 . 4687389E-04 

6 

1 . 9249510E+07 

1 . 1354126E+17 

-6 . 1246793E+05 

7 

- 

-5 . 8887562E+18 

1 . 8533381E+07 

8 

- 

1 . 2840182E+20 

-2 . 7071295E+08 

9 

- 

- 

1 . 9286742E+09 

10 

- 

- 

-5 . 3094692E+09 
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A.  4  Green's  Functions  Describing  Fields  Generated  by  Sources  on 
Interface  3b 

A. 4.1  Green's  Functions  for  Fields  on  Interface  2b. 


Cq32<*)  *  Gq23<*) 

The  coefficients  for  GbA32(R)  are  the  same  as  those  listed  in 
Table  A-9 . 

The  coefficients  for  Gbq32(R)  are  the  same  as  those  listed  in 
Table  A-10. 
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A. A. 2  Green's  Functions  for  Fields  on  Interface  3b. 


■'A33  ' 
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A*2b 

r  h‘2)(ai?)a 
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Jc  -Dfl —  A33 

nbxx 
~  °A33 

(R) 

1 

r  Hq2*  (\R)  X 

where 


^A33  “  /Jbl2u2bcos^1[LI2b(  ^2b  _^lb  )  ] 

+  ulbcoth(blbulb)  sinh[u2b(b2b-blb)] 

^q33  “  {u3b  +  cos^2[u2b  (^2b  _^lb  )  ]  ^bl3ulbtan^(^lbulb)  }  u 


2 

2b 


(ulb  +  Mbl2ebl2ulb)u3b 
*  M b23  €bl2ulbu2bcot^  (blbulb) 


jsinh2[u2b  (b2b-blb)  ] 


r 


2  2 
^b23ulb  +  A*bl3€bl2U2b 


ebl2COth(^lbulb) 
+  /lb12tanh(blbulb) 


fulbu3b 


.  u2b 

~T 


iinh[2u2b(b2b-blb)] 


(A— 20) 


(A-21) 


A-16 


GjSaW 

where 

rA3  "  y 


<33  (R) 
where 

Tg3  "  y 


A  A.  2.1  Asymptotic  Forms  for  R  -*  0. 


^2b  /  f  *c 


<33 *)  "  tA3 


|*RAe  J0(X) 

Jo  E 


dx  + 


(A-22) 


1  ♦  Mbl2 
MbX2  +  Pb23  +  /^b!3 


(A-23) 


j  W  2  £  3b 


{j;c 


0(x) 

IT" 


dx 


Tq3 

~F 


} 


(A-24) 


1  *  £b!2 
£b!2  +  £b23  +  f  b!3 


(A-23) 


and  the  integrals  IbA33(AR)  and  Ibq33(AR)  are  the  integrands  of  the 
original  Sommerfeld  integrals  with  H0<Z)(AR)  replaced  by  J0(AR) . 


Table  A-ll  Coefficients  for  GbA33  Polynomial  Model 


Coefficient  # 

Real 

Imaginary  #1 

Imaginary  #2 

0 

2 . 9895392E-04 

2 . 6798747E-10 

-3 . 2193755E-07 

1 

-3 . 9713808E-01 

-7 . 5940224E-09 

1 . 1728108E-07 

2 

2 . 7003667E+02 

1 . 8608318E-07 

1 . 7359691E-04 

3 

-1.0882422E+05 

1 . 5038343E-14 

1.4780323E-03 

4 

-4 . 2349631E+09 

-3 . 2242885E-04 

-9 . 3812439E-02 

5 

3 . 9979300E+11 

-3 . 1416726E-11 

9 . 3841354E-01 

6 

-2.0878692E+13 

3.4534456E-01 

-2 . 9073252E+00 

7 

2 . 7177121E+07 

-8 . 3578282E+00 

- 

8 

4 . 6215527E+14 

8 . 8658536E+01 

- 

9 

- 

-4 . 6285555E+02 

- 

10 

- 

9 . 7745550E+02 

- 

Table  A-12  Coefficients  for  Gbq33  Polynomial  Model 


Coefficient  # 

Real 

Imaginary  #1 

Imaginary  #2 

0 

9 . 6785892E-01 

-7 . 1240104E+02 

-7 . 3279954E-04 

1 

3 . 2169259E-01 

9 . 5940154E+03 

8 . 4236082E-02 

2 

-1.2286658E+03 

-6 . 5262043E+08 

-3 . 2371602E+00 

3 

-6 . 8608196E+03 

2 . 6293123E+11 

3 . 4357684E-07 

4 

8 . 4780207E+05 

-6 . 5657068E+13 

9 . 1330452E+02 

5 

-1 . 0443582E+07 

1 . 0231241E+16 

—5 . 4107447E-04 

6 

3 . 7986523E+07 

-9 . 6588315E+17 

-2 . 6362783E+05 

7 

- 

5 . 0443342E+19 

9 . 1863973E+06 

8 

- 

-1 . 1166002E+21 

-1 . 8340614E+08 

9 

- 

- 

1 . 6197594E+09 

10 

- 

- 

-5 . 0799782E+09 
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Appendix  B.  Evaluation  of  the  Analytic  Asymptotic  Green's  Function  Term 


This  Appendix  presents  details  of  evaluating  the  integrals  of 
equation  (3-36)  which  is 


HR) 


dsf  dx 


(3-36) 


where  the  contour  Cxj  is  defined  by  x  e  [-a/2,  a/2]  and  the  surface  SXJ 
is  defined  by  x'  e  [-a,  a]  and  y'  €  [-b/2  ,  b/2] .  In  this  instance,  the 
variable  R  -  [ (x  -  x')2  +  (y^  -  y')2]'.  Over  a  y-directed  current  cell 
the  integration  is  defined  as  follows: 


™  •  1 1 

where  CyJ  is  defined  by  y  t  [-b/2,  b/2]  and  the  surface  SyJ  is  defined  by 
x'  t  [-a/2,  a/2]  and  y'  t  [-b,  b]  and  R  -  [  (Xj  -  x')2  +  (y  -  y')2]*. 


i  U  _  i  y'  i 

T?  - 5 - 


ds'  dy 


(B-l) 


B-l 


Equations  (3-36)  and  (B-l)  are  identical  with  the  exchange  of 
variables  a  and  b  except  for  the  constants  Xj  and  y^  which  will  drop  out 
of  the  equations  after  an  exchange  of  variables.  Since  the  expressions 
over  both  x-directed  and  y-directed  cells  are  equivalent,  only  the 
integration  of  (3-36)  will  be  detailed. 

The  first  step  of  the  integration  is  to  complete  the  integration 
with  respect  to  y'  ,  to  do  so  let  u  ->  ym  -  y'  therefore  du  -  -dy'  ,  and 
equation  (3-36)  then  becomes 


HR)  -  2 


1 

y  (x-x')2  *  u2 


a 


du  dx'  dx 


(B— 2) 


where  the  fact  that  the  integrand  is  even  with  respect  to  u  has  been 
used  to  half  the  integration  interval  and  double  the  result. 
Completing  the  integration  with  respect  to  u  produces 


UR)  -  2 


- lb/2 

ln[u  +  \j(x-x'  )2  +  u  J0  dx' dx 


(B-3) 


After  evaluating  the  integrand  for  the  limits  of  u,  the  expression 
for  I(R)  can  be  written 


HR)  -  2 


JL 


i .  -*j  1 


f  (x)  dx' 


(B-4) 


B-2 


where  the  function  f(x)  is  defined  by  the  integral 


*r«/2 

f(x)  -  f  In 


-  J 

-«/2 


1*, 

(x-x7)2  +  -k 

(x-x7  )  1 

dx 


(B-5) 


To  alleviate  the  difficulties  of  the  absolute  value  argument  in 
the  denominator  of  equation  (B-5) ,  the  interval  of  integration  can  be 
broken  into  two  separate  sections  such  that  x'  <  x,  and  x'  >  x. 
Therefore  the  integral  of  equation  (B-5)  can  be  written  as  follows: 


fix)  - 


t! 

i 


In 


xr-«/2 


b 

> 

(x-x')2  +  ^ 

(x7  -x) 

Xj*a/2 


I 


In 


VJ 

(x-x7)2  + 

(x-x7  ) 

dx 


dx 


(B-6) 


Let  the  first  integral  of  (B-6)  be  designated  fx(x)  and  the  second 
as  f2(x) .  To  evaluate  these  integrals,  the  change  of  variables 
(x'  -  x)  -  bu/2  is  used  in  f:(x),  and  (x  —  x')  -  bu/2  is  used  in  f2(x) . 
With  this  substitution,  the  integrands  become  the  logarithmic  identity 
for  the  inverse  hyperbolic  cosecant  of  u,  and  the  integrals  have  the 
form 


B-3 


( B— 7 ) 


f,(x)  -  ^  J*  csch'1 (u)  du 

where 


2 

3 

2 

3 


s  -  1 

s  -  2 


(B-8) 


The  result  for  the  integration  of  equation  (B-7)  is  obtained  from 
an  integration  table  [7:375],  and  then  evaluated  for  the  interval  of 
integration  producing  the  results 


fi(x)  -  [  Aj  csch'1(A1)  +  sinh‘1(A1)] 

f2(x)  -  [ A2csch'l(A2)  +  sinh'1(A2)| 


(B-9) 


Substituting  equation  (B-9)  into  (B-4)  the  expression  for  I(R)  is 
now  given  by 


XJ’« 

I  X1  -X,  I 

A:  csch'1(A1)  +  sinh*1(A1) 

■< 

I(R)  -  b> 

f 

1  *  1  J 1 

dx* 

xrm 

V 

a 

♦  A2  csch'1(X2)  *  sinh*1^) 

* 

(B-10) 
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Now,  consider  the  terms  with  arguments  of  Ax,  and  define  the 
integral  with  the  csch_1(A1)  term  as  I1(x'),  the  integral  with  the 
sinh_l(A1)  as  I2(x'),  the  csch_1(A2)  term  as  I3(x'),  and  the  integral  with 
the  sinh_1(A2)  as  I*(x').  Let  u  -  x'-  x3 ,  then  the  components  of  I(R) 
become 


a 

Ij(u)  -  4j 

1-Jiiil 

L  a  J 

r 

I2(u)  -  2b 

•  1  _  1  U  1 

L  a 

a 

I3(u)  -  A  f 

l  -  JiLl 

4 

a 

r 

I,(U)  -  2b 

1  -  1  u  1 
a 

( B— 1 1 ) 


where  the  fact  that  all  four  integrals  are  even  with  respect  to  u  has 
been  used  to  half  the  interval  of  integration  and  double  the  result.  In 
the  first  and  third  integral,  the  2/b  factor  has  been  combined  with  the 
b  outside  the  integral  in  equation  (B-10) .  Now  the  original  integral  is 
given  by  I(R)  -  AIj(u)  +  2bl2(u)  +  4I3(u)  +  2bl4(u). 

To  find  I1<u),  let  v  -  2(u  +  a/ 2)/b,  then  I2(u)  -  Ia(v)  where 


Ii(v)  -  ^j2  JT  -  Y^jvcsch'HvJdv 


B-5 


(B-12 ) 


Breaking  I^v)  into  two  parts,  the  integration  results  for  the 
parts  are  found  in  the  integral  tables  [7:375],  and  Ij(v)  can  be  written 
as  follows: 


I i(v)  - 


2\v2  csch‘1(v)  sinh'1(v)  ] 


■\v-3r 


b 

Ti 


csch_1(v) 


+  \  sinh'^v)  -  — -v  = 
3  3  fizTT 


(B-13) 


where 

r  -  a/5 

and  the  absolute  value  signs  in  the  integration  formulas  from  the  tables 
can  be  ignored  since  v  is  always  positive  for  v  e  [r,  3r ] . 

For  I2(u) ,  use  the  same  change  of  variables  u  to  v  and  obtain 


I2(v)  - 


fb)  V 

f  3  _  V  ' 

hi  J 

h  77  J 

sinh'^v)  dv 


(B-14) 


I2(v)  can  be  separated  into  two  integrals  and  the  formulas  from 
the  integral  table  [7:374]  give 


For  I3(u)  and  I*(u),  the  substitution  v  -  2(a/2  -  u)/b  is  used  and 
the  interval  of  integration  is  given  by  v  t  [-r,  r].  Since  v  can  be 
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v-3f 


b 

Ti 


Jv2  *  1 


V  +  1 


value  of  v  must  be  retained  in  the  integral  formulas  and  the  results  are 


I3(v) 


-w 


-VV-T 


3 

% 


TS 


v2  csch_1(v)  +  -pi-p  sinh_1(v) 


^ csch-Hv ) 


TT7T 


sinh'Mv)  - 


sjv2  + 1 


(B-16) 


I*(v) 


-i  v  sinh‘1(v)  -  ^v2  +  1 


b 

Ti 


^  jsinh'Mv)  -  ^ 


(B-17) 


Evaluating  equations  (B-14)  through  (B-17)  at  the  defined  limits, 
combining  the  results,  and  using  the  identities  sinh_1(-x)  -  -sinh_1(x) 
and  csch_1(-x)  -  -csch_1(x)  the  expression  given  in  equation  (3-37)  is 
obtained . 
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